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ABSTRACT 

In  Project  Early  Rise,  the  University  of  Alberta 
recorded  arrivals  at  distances  of  1100  to  2000  kilometers 
from  the  source.  At  some  locations  a  cluster  of  one  vertical 
and  two  horizontal  components  of  ground  motion  were  recorded. 

A  very  fast  and  efficient  program  for  Butterworth 
band-pass  filter  has  been  written  and  was  used  to  enhance  the 
data  without  introducing  any  phase  distortion.  This  algorithm 
based  on  the  Butterworth  filter  has  given  very  good  results 
in  various  geophysical  problems  involving  time  sequence 
analysis . 

Theoretical  spectral  ratios  have  been  compared  with 
the  spectral  ratios  of  field  records  from  Project  Early  Rise. 
The  theoretical  spectral  ratios  for  a  southern  Alberta  crustal 
model  with  five  layers  has  been  calculated  using  the  Thomson- 
Haskell  matrix  formulation.  The  basic  Alberta  model  has  been 
modified  by  changing  the  thickness  of  various  layers  to  ob¬ 
tain  a  better  fit  between  theoretical  and  experimental 
spectral  ratios  at  various  locations  where  data  was  available. 
From  this  limited  amount  of  comparisons  it  is  found  that  the 
spectral  ratio  favors  a  crustal  thickness  of  about  43  kilo¬ 
meters  in  Saskatchewan  and  increasing  to  48  kilometers  on  the 
Alberta-Saskatchewan  boundary  to  a  possible  maximum  of  53 
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kilometers  near  Viking,  Alberta.  No  uniqueness  can  be  claimed 
for  these  models  and  it  is  probable  that  further  calculations 
would  yield  models  which  give  a  better  fit  over  the  frequency 
range  that  was  recorded. 

Comparison  of  field  results  were  also  made  using 
theoretical  seismograms  which  were  calculated  by  a  fast  Fourier 
program  which  inverted  the  transfer  function  obtained  with 
a  Thomson-Haskell  matrix  formulation.  It  was  noticed  that  at 
certain  locations  the  theoretical  seismograms  agree  with  the 
field  records  while  at  other  locations  the  comparison  was  not 
good.  Possibly,  the  crustal  model  will  need  major  modifica¬ 
tion  at  these  locations  or  the  assumption  of  plane  layering 
In  the  model  calculations  is  not  justified. 
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CHAPTER  1 

SEISMIC  INVESTIGATION  OF  THE  CRUST 

1.1  Introduction 

An  investigation  of  the  deep  structure  of  the 
Earth  will  help  us  to  interpret  its  historical  evolution, 
the  nature  of  some  destructive  events,  and  it  will  even  give 
us  a  better  understanding  about  other  planets.  In  the  study 
of  the  Earth’s  interior  refraction  and  reflection  seismology 
have  been  two  of  the  most  powerful  tools.  From  the  varia¬ 
tion  of  velocity  with  depth  the  Earth  has  been  divided  into 
a  number  of  zones  which  are  separated  by  first  or  second 
order  discontinuities.  The  uppermost  part  of  these  zones 
is  called  the  crust  and  it  has  an  average  thickness  of  35  km 
under  the  continents  and  5  km  under  the  oceans.  The 
Mohorovicic  (M)  discontinuity  separates  the  crust  from 
another  zone  called  the  mantle  having  different  elastic  para¬ 
meters  . 


The  crust  has  been  found  to  consist  of  many  dis¬ 
crete  layers  whose  properties  change  from  place  to  place. 
This  layered  medium  produces  a  reverberation  which  makes  it 
difficult  to  interpret  deeper  arrivals  from  the  mantle. 

This  aspect  of  the  problem  will  be  investigated  at  length  in 
this  thesis.  Apart  from  this  effect  a  study  of  a  layered 
system  as  a  linear  filter  is  useful  because  of  the  economic 
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importance  of  the  crust. 

1.2  Previous  Work. 

Observations  of  the  effects  of  the  crust  on 
seismic  reverberations  were  begun  by  Imamura  (1929)  who 
studied  terrestrial  vibrations  induced  by  the  arrival  of 
seismic  waves  from  earthquakes.  He  observed  that  the 
period  of  the  waves  varied  with  different  soils  and  geolo¬ 
gical  strata.  Nasu  (1931)  continued  studying  vibrations 
of  the  Earth  and  concluded  that  earthquake  motion  at  short 
periods  (less  than  1  sec.)  is  much  influenced  by  the 
nature  of  the  soil  and  he  explained  the  oscillations  as  due 
to  a  superficial  layer  being  excited  by  a  sudden  arrival 
of  an  impulse  originating  from  an  earthquake.  Suzuki  (1932) 
observed  changes  in  the  period  of  waves  with  the  angle  of 
incidence.  His  observational  study  showed  that  beyond  cer¬ 
tain  epicentral  distances  the  angle  of  incidence  does  not 
change  with  distance.  Gutenberg  (1957)  studied  earthquake 
damage  at  different  locations  as  a  function  of  the  elastic 
properties  of  the  soil.  Nuttli  (1964)  observed  that  for 
long-period  P  waves  generated  by  earthquakes  of  magnitude 
6,  the  amplitude  of  ground  motion  was  not  effected  by 
crustal  layering. 

The  theoretical  response  of  such  a  layered  system 
has  been  calculated  by  different  authors  under  certain 
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simplifying  assumptions.  For  example,  Sezawa  (1930)  cal¬ 
culated  the  free  oscillations  of  a  single  layer  excited  by 
a  dilatational  pulse  of  a  purely  plane  type  propagated  ver¬ 
tically  upwards  in  an  elastic  solid  medium.  Later  Sezawa 
and  Kanai  (1932  a,b,  1937)  studied  transmission  of  seismic 
waves  of  oscillatory  type  and  of  finite  extent  through  a 
strata  for  obliquely  incident  waves.  Theoretical  seismo¬ 
grams  were  compared  by  changing  the  thickness  of  the  layers 
and  the  frequency  of  the  waves.  Kanai  (1952,  1953,  a,b) 
in  a  series  of  studies,  calculated  the  oscillations  of 
doubly  stratified  visco-elastic  layers  excited  by  seismic 
waves  at  vertical  incidence  and  the  relation  between  the 
amplitude  at  the  surface,  the  impedance  ratio  of  the  two 
media,  the  coefficients  of  viscosity  and  the  thickness  of 
the  surface  layer.  He  showed  that  as  the  velocity  of  a 
weak  layer  decreases,  the  amplitude  of  the  seismic  motion 
at  the  free  surface  increases. 

Thomson  (1950),  Haskel  (1953)  and  Dorman  (1962) 
generalized  the  theory  to  yield  the  response  of  a  layered 
system  to  plane  waves  at  any  angle  of  incidence.  They 
expressed  the  solution  in  terms  of  products  of  four  by 
four  matrices  whose  members  are  a  function  of  the  parameters 
of  each  layer;  the  angle  of  incidence  and  the  wavelength 
of  the  incident  harmonic  wave  .  Thomson  (1950)  made  the 
original  formulation  of  propagation  of  plane  waves  through 
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a  system  of  plane  parallel  layers.  Haskell  (1953)  corrected 
an  error  in  Thomson's  work  and  discussed  the  application 
of  the  method  to  the  dispersion  of  Love  and  Rayleigh  waves. 

In  1962  Dorman  formulated  the  solution  for  the  boundary  con¬ 
ditions  at  a  boundary  between  a  liquid  and  a  solid.  Later 
Haskell  (1962)  gave  the  transfer  function  for  the  vertical 
and  horizontal  components  of  ground  motion. 

Fast  digital  computers  have  made  it  practical  to 
use  the  Thomson-Haskell  matrix  formulation  for  calculating 
the  response  of  a  media  with  any  number  of  layers.  Phinney 
(1964)  compared  the  ratio  of  long-period  spectra  of  vertical 
and  horizontal  seismograms  for  large  earthquakes  with  the 
ratio  of  the  vertical  to  horizontal  transfer  function,  using 
the  Thomson-Haskell  matrix  formulation.  Hanon  (1964)  first 
attempted  the  calculation  of  a  synthetic  seismogram  by 
taking  the  inverse  Fourier  transform  of  the  transfer  function 
for  different  models.  Fernandez  (1965)  calculated  spectral 
ratios  for  long-period  data  giving  model  curves  to  determine 
crustal  structure.  Recently  Leblanc  (1966)  and  McCamy 
(1967)  applied  the  transfer  function  and  spectral  ratios  to 
short  period  data  in  an  attempt  to  obtain  crustal  structure. 

In  this  study,  the  transfer  function  for  an 
ALBERTA  crustal  model  has  been  calculated  using  the  Thomson- 
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cal  to  horizontal  transfer  function  was  compared  with  the 
ratio  of  the  spectra  of  vertical  and  horizontal  seismograms. 
A  time  synthesis  of  ground  motion  was  obtained  for  a 
Sin  X/X  shaped  pulse.  Comparisons  were  carried  in  the 
time  domain. 

1.3  Experimental  Data. 

Project  Early  Rise  was  planned  to  study  the  upper 
mantle  under  North  America  by  seismic  waves  generated  by 
chemical  explosives  detonated  in  Lake  Superior.  The  Univer¬ 
sity  of  Alberta  recorded  the  arrivals  at  a  distance  of  1100 
to  2000  kilometers  from  the  source.  Arrivals  are  recorded 
also  by  a  400  kilometers  arc  line  from  Cold  Lake  to  Bow 
Island,  Alberta.  Nine  stations  on  this  arc  line  have  a  dis¬ 
tance  of  1675  kilometers  from  the  source.  The  United  States 
Geological  Survey  extended  this  profile  from  350  to  1100 
and  2000  to  3500  kilometers  from  the  shot  to  the  Alaska- 
Yukon  border.  The  seventy-four  locations  are  labelled  from 
WC  1  to  WC  74  (Figure  1.1). 

In  this  project  the  University  of  Alberta  had 
two  mobile  recording  units  equipped  with  tape  recorders  and 
two  stationary  observatory  sites  also  recording  on  magnetic 
tape.  Each  mobile  unit  had  9  Texas  instrument  S-36 
2  c.p.s.  seismometers,  12  channel  V.L.F.  amplifiers,  an 
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Figure  1.1.  Project  Early  Rise 
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RS-8U  24-trace  oscillograph  with  200  c.p.s.  galvanometers 
and  a  precision  instrument  model  PS  207  F.M.  magnetic  tape 
recorder.  Occasionally  a  cluster  of  one  vertical  and  two 
horizontal  components  of  ground  motion  was  recorded  using 
Hall  Sears  (HS-10)  2  c.p.s.  seismometers. 

Most  of  the  recordings  were  made  with  high-cut 
filters  on  8  c.p.s.  The  low  frequency  cut-off  is  not  adjus¬ 
table  but  falls  off  rapidly  below  2  c.p.s.  Attenuation  of 
arrivals  was  usually  set  to  18  db  on  the  amplifier  panel. 

A  chronometer,  WWV  receiver  and  standard  broadcast  receiver 
were  used  for  timing  purposes. 

The  accuracy  of  the  results  obtained  about  the 
Earth’s  interior  depend  upon  the  reliability  of  the  data. 

The  use  of  precisely  located  and  timed  explosion  is  preferred 
to  studies  employing  earthquake  generated  waves.  Current 
recording  systems  using  magnetic  tape  allows  us  to  enhance 
the  data  with  various  data  processing  techniques  on  the 
digital  computer.  Useful  data  was  obtained  during  Project 
Early  Rise  and  in  Chapter  2  attempts  have  been  made  at 
enhancing  this  data  with  distortionless  digital  band-pass 
filters.  The  results  were  then  used  in  Chapter  3,  to  evalu¬ 
ate  the  usefulness  of  present  day  theoretical  studies  as 
applied  to  the  calculation  of  theoretical  seismograms. 
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CHAPTER  2 

ZERO  PHASE  DIGITAL  FILTERS 

2.1  Introduction. 

Zero  phase  filters  are  most  important  because  they 
are  the  best  filters  in  a  least-mean  square  sense  for  ex¬ 
tracting  signals  whose  spectra  overlap  slightly  with  that 
of  noise.  It  is  possible  to  obtain  data  with  no  phase  dis¬ 
tortion  by  digital  filters  since  past  and  future  values  of 
the  data  are  available  in  the  memory  of  the  computer.  This 
is  one  of  the  major  advantages  of  digital  over  analog  filters. 


Figure  2.1. 


Ideal  low-pass  filter 
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An  ideal  filter  is  the  one  with  unit  response  in 
the  desired  region  and  zero  elsewhere  (Figure  2.1)  in  the 
frequency  domain.  An  approximation  to  such  a  filter  can 
be  realized  by  a  number  of  methods  which  are  discussed  by 
Baum  (19^8),  Guillemin  (1957),  Ormsby  (1961),  Weinberg 
(1962),  Kaiser  (1963),  Holtz  and  Leondes  (1966),  Stalling 
(1966)  ,  Shanks  (1967)  and  Wood  (1968). 


is 


In  the  digital  domain  the  unit  impulse  function 


1  t  =  0 
0  t  ?  0 


2.1 


The  output  of  a  filter  whose  input  is  a  unit  impulse  is 
called  the  impulse  response:.  The  filter  operation  is  then 
defined  as  the  convolution  of  the  input  data  with  the  im¬ 
pulse  response: 


n 

y  =  7  f .  x  2.2 

J  n  L  1  n-i 

i=0 

where  (xQ,  x-^  x2,  ...  xm)  are  the  m+1  values  of  the 
input  series; 

(fQ,  f  ,  f2  ...  fR)  are  the  n+1  values  of  impulse 


response ; 


39UlBV  I+m 

I+fl 
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(yQ5  y1 ,  y2  •••  yn+m)  are  n+m+1  values  of  the  out¬ 

put  series . 


For  a  filter  to  be  stable  the  sum  of  the  squares 
of  the  impulse  response  must  be  bounded  i.e., 

00 

i  iftr  < »  2-3 

t  =-00 

A  filter  is  called  realizable  if  f^_  is  one-sided  so 
that  no  response  can  occur  prior  to  its  stimulus.  In 
equation  2.3  the  lower  limit  can  be  set  to  t=0  for  a 
realizable  filter. 

For  digital  computation  and  analysis  of  filter 
stability  a  z-transform  is  used  (Jury,  1964),  since  con¬ 
volution  in  the  digital  time  domain  is  equivalent  to  multi¬ 
plication  in  the  domain  of  the  z-transform.  The  z-transform 
of  a  continuous  function  f(t)  which  has  a  sampling  inter¬ 
val  of  At  is 

00 

F(z)  =  l  f(nAt)  Zn  2.4 

n=-°° 

where  n=0,l,2  ...  and  z  is  a  complex  variable.  The 
z-transform  of  a  realizable  stable  system  can  be  written 
as 


00 

m9oov,o  jl&scta  @ Lcffis.'-Le^'i  fi.  o 


00 
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Y(z)  =  l  y ( t )  zt 
t  =  0 


2.5 


where  y(t)  is  the  impulse  response  of  the  filter.  The 

—  A  t  s 

z-transform  can  also  be  obtained  by  replacing  e  or 

e  with  z  in  the  Laplace  and  Fourier  transform 

(Jury,  1964).  Therefore  z-transforms  can  also  be  expressed 
as 


z 


-At  s 
e 


-icoAt 

z  =  e 


2 . 6 

2.7 


Here  z  represents  the  delay  of  a  signal  one  sample  inter¬ 
val  and  z2  represents  a  delay  of  two  sample  intervals 
(Robinson,  1964),  etc.  Using  the  equation  2.7  in  the  equa¬ 
tion  2.5  for  At  =  l,  that  is,  on  the  unit  circle  in  the  z-plane, 
it  can  be  written 

00 

Y(oo)  =  £  y(t)  e  -7T  <  (jO  <  TT  2.8 

t=0 

The  function  Y(oo)  is  called  the  transfer 
function  of  the  filter.  In  polar  form  it  can  be  written 


as 


(  ;  j  dO  9.  C  '.3  '  ' >  ■ :  ~ “ 


lo  ’  '9... i  s  so  e  .  ri 


V  .  ;  i  "  'Jf  3  £. 
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Y(co)  -  |  Y ( ca )  |  e10(“) 


2.9 


where  |  Y ( a) )  |  is  called  the  Gain  and  0(co)  the  phase 
shift.  The  delay  of  the  system  is  defined  as 


dQ(co) 

dm 


2.10 


2.2  Recursion  filters. 

The  filter  operation  in  the  z-domain  can  be 

written  as 


Y ( z  )  =  F(z) •  X(z) 


2.11 


where 


X(z)  =  x0  +  x^z  +  x^z2  + 


+  X 


m 


z 


m 


F(z)  =  fQ  +  fxz  +  f2z2  +  .  +  Vn 

Y(z)  =  y0  +  yxz  +  y2z2  +  .  +  ym+nzm+n 


2.12 

Often  the  transfer  function  of  a  filter  can  be 
expressed  as 


F(z) 


A(  z) 

B  ( z ) 


an  +  an z  +  a0z2  +  ....  +  a  zn 
0_ 1_ 2_ n 

o  m 

bn  +  bnz  +  bnz2  +  ....  +  b  z 
0  12  m 

2.13 
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Such  a  transfer  function  can  be  programed  more  efficiently 
by  a  recursion  relation  involving  a  feedback  loop.  As  an 
example  consider  an  impulse  response  with  the  form 


F(z) 


aQ  +  a±z 


1  +  b-^z  +  b^z 


2.14 


From  equation  2.11  the  output  will  be 


Y(z)  =  F(z) •  X ( z ) 


aQ  +  a]_z 


1  +  b,  z  +  b~z 


X  ( z ) 


1 


2.15 


Multiplying  both  sides  of  equation  2.15  by  the  denomina¬ 
tor  and  arranging  terms  we  obtain 


Y(z)  =  [aQ  +  a^l-XCz)  -  z  Y(z)[b1  +  b2z]  2 . 16 


Equation  2.16  can  be  evaluated  by  following  the  steps 
indicated  in  the  flow  diagram  of  figure  2,2. 


Figure  2.2.  Evaluation  of  the  recursive  equation  2.16. 
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A  general  recursion  equation  representing  equation  2.13 
can  be  written  as 


14 


n 


m 


y 


n 


y  a .  x 

L  1  n-i 


2.17 


i=0 


where  the  coefficient  bQ  is  assumed  to  be  1.0  for 
convenience  in  equation  2.13.  The  term,  bQ  can  always  be 
made  equal  to  1.0  by  division.  An  algorithm  based  on 
the  recursive  relation  in  equation  2.17  can  save  a  great 
deal  of  computer  time  and  is  more  precise  than  a  standard 
convolution  (equation  2.2).  Equation  2.2  requires  many 
more  multiplications  and  summations  and  requires  a  large 
number  of  filter  coefficients,  f.  . 

2.2.1  Synthesis  of  recursion  filters  in  the  z-plane. 

It  is  useful  to  analyze  rational  filters  in 
terms  of  roots  of  the  z-transform  of  the  numerator  and 
the  denominator.  The  roots  of  the  numerator  are  called 
zeros  of  a  filter  while  the  roots  of  the  denominator  are 
called  poles  of  a  filter.  A  plot  of  the  location  of  the 
poles  and  zeros  in  the  complex  z-plane  may  be  used  to 
predict  the  stability  of  a  filter  operator  (Figure  2.3). 
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Figure  2.3.  z-plane  representation  of  poles  (x)  and 

zeros  (o)  . 

The  amplitude  and  phase  response  of  a  digital 

filter  can  be  found  by  evaluating  on  the  unit  circle  |z|  =  1. 

From  equation  2.7  the  point  z  =  1.0  +  i  0.0  corresponds 

to  zero  frequency  and  z  =  -1.0  +  i  0.0  corresponds  to 

f  =00  / 2tt  =  — —  s  the  Nyquist,  folding  or  alias  frequency. 
n  n  2At 

Frequencies  are  linearly  distributed  around  the  unit  circle 
between  zero  and  the  Nyquist  frequency  (Figure  2.3). 
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The  z-plane  can  be  used  to  design  simple  filters  by 
locating  poles  and  zeros  which  are  known  (Truxal,  1955). 

The  poles  must  be  located  outside  the  unit  circle  since  a 
pole  inside  the  unit  circle  causes  instability  (Treitel 
and  Robinson,  1964).  The  amplitude  response  of  a  filter 
at  any  desired  frequency  is  related  to  the  ratio  of  the 
distance  of  the  zeros  to  the  distance  of  the  poles  from  the 
frequency  point  on  the  unit  circle  in  the  z-plane  (Figure 
2.3)  . 


2.3  The  Bilinear  z-transform. 


Golden  and  Kaiser  (1964)  have 
algebraic  transformation  which  converts 
fer  function  into  one  which  can  be  used 
It  eliminates  aliasing  errors  which  are 
standard  z-transform  (equation  2.6). 
z-transformation  is  defined  as 


employed  an 
a  continuous  trans- 
on  sampled  data, 
present  in  the 
The  bilinear 


2  ,  ,  / s At  s 

at  tanh  (_2_) 


2.18 


using  equation  2.6  we  have 

_  2  H(eiss(At)  -  eiss(At)) 

sd  At  j^g^sCAt)  +  e*5S(At) ) 


2.19 


or 
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s  =  1-z 

d  At  1+z 


2.20 


Non-linear  warping  caused  by  this  transformation 
can  be  seen  by  setting  sd  =  io3  in  equation  2.18  as 
illustrated  in  figure  2.4. 


t  _  2  ,  (wAt) 

wd  At  tan  2 


2.21 


Since 


tanh  x  =  -i  tan  ix 


—  TT  ^  7T 

——  <  03  <  - - 

At  At 


where 


ood  =  Deformed  angular  frequency  used  to  calculate 
s^  and  the  bilinear  z-transform. 

o)  =  Angular  frequency  in  the  original  transfer 
function. 

A  correction  must  be  made  by  equation  2.21  before  the 
poles  and  zeros  of  a  filter  are  determined. 


5  1  f  r  At.  r  •'  .  -c  ‘  -  a  ‘  ^  Y.c  *r-  1  - 


£i  1  ,  c  i  _L_  =  t. 
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Figure  2.4. 


Non-linear  warping  of  frequency  scale  in  the 
bilinear  z-transformation . 


This  transformation  maps  the  entire  oo^  into  the  horizon¬ 
tal  strip  bounded  by  Nyquist  frequencies: 

(jo-,  =  -ioo  /2  and  oo0  =  +ioo  /2  2.22 

Is  2  s 


The  bilinear  transformation  is  similar  to  the  standard 
z-transform  at  low  frequencies  as  shown  in  figure  2.4. 
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2.4.  Butterworth  filters. 

To  approximate  an  ideal  filter,  it  is  necessary 
to  truncate  the  impulse  response  before  it  becomes  negli¬ 
gibly  small,  so  that  computation  time  is  not  excessive. 
This  causes  a  ripple  or  oscillation  called  Gibb’s  effect 
which  is  a  very  serious  defect  in  these  types  of  filters. 
To  overcome  this  difficulty  a  Butterworth  function  is  used 
because  it  has  a  flat  response  over  the  pass  region  where¬ 
as  other  possible  functions  have  ripples.  It  can  be  shown 
that  the  Butterworth  function  is  the  best  approximation 
to  an  ideal  filter  (Baum,  1948). 

The  Butterworth  function  (Butterworth,  1930)  is 

given  by 

|Y(S)  | 2  =  - - -  2.23 

1  +  ?52n 

The  symbol  ft  is  used  to  represent  the  normalized  low- 
pass  frequency.  For  a  cut-off  frequency  ooL  ,  the  norma¬ 
lized  frequency  will  be 


ft 


oo 


2 . 24 


The  equation  2.23  is  maximally  flat  over  the  pass  region 
and  a  larger  value  of  n  gives  a  greater  rate  of 
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attenuation  (Figure  2.5). 


Figure  2.5.  Square  of  the  transfer  function  for  a  norma¬ 
lized  low-pass  filter  using  a  Butterworth 


function . 
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2.4.1.  Poles  and  zeros  of  a  low-pass  Butterworth  filter. 

The  poles  and  zeros  are  obtained  from  equation 
2.23  by  making  the  frequency  complex  as  p  =  p  +  ifi  with 

p  =  0  : 


YT  (p) 


1  =  1 _ 

1  +  (|)2n  1  +  (-nVn 


2.25 


Since  the  numerator  is  constant  there  are  no  zeros  in  a 
low-pass  filter  except  at  infinity.  The  poles  occur  at 
places  where  the  denominator  is  zero: 


1  +  (-l)n  p2n  =  0 


2.26 


or  if  we  multiply  both  sides  by 


(-l)n 


2n  , 

P  + 


(-Dn  =  0 


2.27 


The  poles  will  be  arranged  around  the  unit  circle  in  the 
complex  frequency  domain.  Only  the  poles  on  the  left- 
hand  side  of  the  complex  plane  are  taken  to  insure  that 
we  have  stability  and  minimum  phase  (Weinberg,  1962). 

Equation  2.23  can  be  written  as 


1 


Bn(p)  '  Bn(p) 


yl(p) 


2 


2 . 28 
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where  B  (p)  is  called  the  Butterworth  polynomial.  The 
roots  of  the  transfer  function  which  are  on  the  left  side  of 
the  p-plane  are  given  by 


p2v+l 


=  e 


i7T/2n  •  ( 2v+l+n) 


v=0 , 1 , . . . . n-1  2.29 


or  by 


p0„±1  =  -  Sin  ■—(  2  v+1 )  +i  Cos  ~(2v+l)  2.30 


2v+l 


2n 


2n 


The  Butterworth  polynomial  can  be  written  as 


„  /  s  n  1  /  i7T/2n.(  2v+l+n)  s 

B  (p)  =  tt  (p  -  e  ) 

v  =  0 


2.31 


This  can  be  written  (Weinberg,  1962)  in  the  form  of 


n 


Vp)  =  ,  I  ak  p 

k=0 


k 


2.32 


where 


p  =  n  Cos (y-1) y 
k  y=i  Sln  VY 


2.33 


and 


Y  = 


TT 

2n 


a0 


1 
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B 
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For  example  the  poles  of  a  normalized  fourth  order 
low-pass  filter  as  obtained  from  the  equations  above  are: 


P1  4  =  -0.3826834  ± 

i  0.9238795 

2.34 

P2]3  =  -0.9238795  ± 

i  0. 3826834 

2.35 

As  required  by  theory  these  occur  on  the  unit  circle  in  the 
p-plane  (Figure  2.6).  If  the  frequency  is  not  normalized 
the  poles  will  be  outside  the  unit  circle.  For  n=4 
equation  2.28  can  be  written  as 


yl(p)  (p-p1) (p-p2) (p-p3) (P-P4T  2,36 

where  p P  q  n  are  the  four  poles  of  the  Butterworth 
function  for  n=4  .  . 


Figure  2.6.  Poles  of  a  fourth  order  low-pass  Butterworth 


filter . 
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For  a  desired  cut-off  frequency,  oo  ,  the  transfer  func- 

0 

tion  can  be  obtained  by  replacing  p  by 


Pn  =  P  •  “c  2-37 

To  obtain  a  recursive  relation,  it  is  convenient 
to  express  the  denominator  of  equation  2.36  in  terms  of 
second  order  equations: 

YT(p)  =  - i -  2.38 

(p2+bp+c) (p2+dp+e) 

Using  the  bilinear  z-transform,  the  transfer  function 
can  be  written  as 


(1+Fl-z  +  F2. z2) (1+F3-Z  +  F4 • z2 ) 

2.39 

where  F1,F2,F3,F4  are  the  filter  coefficients,  which  can 
be  calculated  by 

2 . c  -  8/At2 

FI  =  -  2.40 

4/At2  +  2b/ At 

4/At2  -  2b/At  +  c 


F2 


4/At2  +  2b/At  +  c 


2.41 


- - 


dA\8  -  a.S 
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4-e  -  8/At2 

F3  =  - -  2.42 

4/At2  +  2d/At 

4/At2  -  2d/ At  +e 

F4  =  -  2.43 

4/At2  +  2d/At  +e 

Equation  2.39  can  be  performed  by  a  recursive  relation  in 
the  computer.  A  fortran  IV  program  has  been  written  for 
IBM  360/67  to  calculate  filter  coefficients  in  equations 
2.40,2.41,2.42,2.43  and  to  process  any  digitized  data  with 
this  filter  (Appendix  2).  As  an  example,  for  cut-off 
frequency  f  =10  c.p.s.  the  p-plane  poles  of  a  low- 
pass  filter  are 


Pl  2  =  -59.4620  ±  i  24.6305 


P 


3,4 


-24.6305  ±  i  59.4620 


2.44 


The  impulse  and  amplitude  response  of  this  filter  is  shown 
in  figure  2.7  and  2.8. 
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Figure  2.7.  Zero  phase  impulse  response  of  a  Butterworth 

low-pass  filter  for  f  =10.0  c.p.s.,  At=l/ll8 

V-* 

and  n=4 


sec  . 
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Figure  2.8.  Zero  phase  amplitude  response  of  Butterworth 

low-pass  filter  for  f  =10.0  c.p.s. 

c 

At=l/ll8  sec.  and  n=4 
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2.4.2  Frequency  transformation  for  high-pass  filter. 

A  high-pass  filter  may  be  constructed  from  a  low- 
pass  filter  by  a  transformation.  Let  the  frequency  in  a 
high-pass  function  be  given  by 


2.45 


s  =  3  +  ico 


where  co  represents  the  normalized  frequency.  It  is 
necessary  to  find  a  frequency  transformation  of  the  type 
p  =  to(s)  which  will  convert  one  type  of  filter  into  another. 
For  a  high-pass  filter  the  transformation  is  found  empiri¬ 
cally  to  be 


1 


2 .46 


by  substituting  p  and  s  we  have 


1 


2.47 


p  +  ift 


3  +  ioo 


or  it  can  be  written  as 


2.48 


If  we  let  3  approach  zero  the  frequency  transformation 


becomes 


.•isdin  sssq-fisxd  iol  noIl£mTL0'’.en£-  v.n  •  ^ 


-wol  £  moil  bedouidanoo  sd  ts^li  £; “  s-- 

b  ni  ^onexps'il  sitf  deJ  .nod  JsiTrcol  insud  £  w  s  >dldl  asf  q 

nnvis  9  nc  tonul  aseq-rl^xri 
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n  =  -  3  2. <49 

03 

t  h 

To  find  the  n  order  high-pass  Butterworth  filter  equa¬ 
tion  2.49  is  substituted  in  equation  2.23: 


Yr(03) 


1 

i  +  (— )2n 

03 


1 


1  + 


2.50 


or 


Yh(u) | 2 


w2n 


1  +  w2n 


2.51 


Figure  2.9* 


Square  of  the  transfer  function  for  an  n^h 
order  Butterworth  high-pass  filter. 
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To  find  the  transfer  function  for  a  high-pass 

filter  with  cut-off  frequency  go  5  the  poles  of  a  low- 

c 

pass  filter  are  converted  to  a  high-pass  poles  using  equa¬ 
tion  2.46: 


2.52 


Equation  2.36  is  used  to  find  the  transfer 

function : 


F  ( s )  =  - — -  2.53 

(s  2+bs+c ) (s  2+ds  +  e) 

Using  the  bilinear  z-transform,  one  can  obtain  the 
impulse  response  in  the  z-plane: 

(1-z)4 

w  =  -  2.54 

n  (1+Fl.z+F2.z2) (1+F3.Z+F4 . z2) 

where  the  filter  coefficients  F1jF2,F35F4  are  given  by 

At  2 • c/2-2 

FI  =  -  2.55 

At/2‘b+At2/4‘C+l 

l-b-At/2+c- At2/4 

2.56 


F2 


b*At/2+c*At2/4+l 


assq-riSlrl  *>  lol  nor'pmJl  *191 11.  ?  •  "  all  f 

-wol  b  lo  39loq  ocU  e 

®  K  I  o'u 


_  ss  W 

<  *s .  t^+s  .  £*+! )  (  s  s  .  S-5+S  .  I^I+I )  . 

aln9XO±llSOO  19  111  J  9/’. 


£-£\o- S1A 
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e  •  At  2/2-2 

F3  =  - - - 

d*At/2+e*At2/4+l 


2.57 


1-d • At/2+e • At  2/4 

F4  =  - - -  2.58 

d • At/2+e • At  2 / 4+1 


where  At  is  the  digitizing  interval. 

A  fortran  IV  program  has  been  written  for  IBM 
process  any  digitized  data  with  the  high-pass 
(Appendix  2).  As  an  example,  the  impulse  and 
response  of  a  high-pass  filter  has  been  shown 
2.10  and  2.11. 


360/67  to 
filter 
amplitude 
in  figure 


2.4.3  Frequency  transformation  for  a  band-pass  filter. 

The  normalized  band-pass  transformation  can  be 
expressed  in  terms  of  a  low-pass  and  high-pass  transforma¬ 
tion  : 

p  =  s  +  i  2.59 

o 

where  p  is  the  complex  frequency  in  the  low-pass  func¬ 
tion  and  s  is  the  complex  frequency  in  the  band-pass 
function. Substituting  complex  s  and  p  values: 


9  ito 


p  =  9  +  ioo  + 


9  2  +  m2 


9  2+oo2 


2.60 
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Figure  2.10.  Zero  phase  impulse  response  of  a  Butterworth 

high-pass  filter  for  f  =10.0  c.p.s., 
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At=l/ll8  sec.  and  n=4 
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Figure  2.11.  Zero  phase  amplitude  response  of  a 

Butterworth  high-pass  filter  for  f  =10.0  c.p.s. 

0 

At=l.ll8  sec.  and  n=4 
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AMPLITUDE  RESPONSE  OF  BUTTERWORTH  HIGH  PASS  FILTER  FOR  10.0  HZ. 
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or 


p  + 


8  (9  2  +go2  +  1  )  _j_ 
8  2+w2 


ioo(a  2  +  go2-1) 

a  2+m2 


2.61 


Letting  the  attenuation  factor  9  approach  zero  the 
frequency  transformation  becomes 


U) 


2.62 


This  is  a  quadratic  equation  having  two  roots,  go  and 

go  : 

2 


GO2  -  HgO  -1  =  0 


2.63 


The  two  roots  are 


j 


2.64 


For  a  unit  cut-off  frequency  in  the  low-pass  domain  i.e., 
U  =  1  ,  the  roots  are 


GO 


1  2 
5 


=  1  +  / 
2  " 


5 


2.65 


or 


Sd.S 


+  =  fit  *  Q 

6  'locfoBl  noxd’Bunbd’ctB  aritf  gnl^sJ 
S9rnoo9cf  noitfBfrnco'lartB'icJ’  yon.  .  •  J3ci'l 

=  a  '  * 


U) 


ts^oon  owcf  gnlvfiri  noi;^BJjp9  oictB^buup  ;  i . 


91B  3  Cf  00*1  ov.xf  9riT 


'ti  :  i  'I:  B8.sq-.JO.  r  l-  ni  v;on  ijjps-r  '  Tio  3  b  )  3 
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a)  =  -0 . 6l8 

1 


2.66 


00  =  1.618 

2 


Setting  n  =  -1  we  obtain  the  conjugate  pair  of  fre¬ 
quencies.  These  four  roots  specify  the  band  width  of  the 
normalized  band  pass  filter  (Figure  2.12). 
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to  a  band-pass  filter. 
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Any  two  arbitrary  frequencies  in  the  band-pass 
domain  may  be  obtained  by  the  appropriate  transformation 
from  a  low-pass  filter.  For  the  normalized  band-pass  we 
have  the  relation 


co  w  =  1  2.67 

1  2 

This  can  be  verified  by  substituting  co  and  co  from 

1  2 

equation  2.64. 

The  normalized  relation  can  be  generalized  to  the 
desired  band-pass  filter  with  unnormalized  upper  and  lower 
frequency  limits  co  and  co  .  The  geometric  mean,  con  , 
which  is  sometimes  called  the  center  frequency  is 

coA  =  /  co  co  2.68 

u  1  2 

or  for  any  cut-off  frequencies  co^  and  co^ 

“0  =  /_vT  2-69 

The  band  width  Aco  is  related  to  the  low-pass  cut-off  as 
seen  in  figure  2.12. 


Aco  =  =  co  -  co 

2 


1 


1 


2.70 


£*-„  s  b  ;  sc:  iis^xlBmnon  Brief  io'? 
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or 

Aw  =  -  0)^  2.71 

The  band  width  will  be  rescaled  by  dividing  by  the  geo¬ 
metric  mean: 


ft 


h 


V“ 


l 


2.72 


Prom  equation  2.69  and  2.72  it  is  possible  to  eliminate 


%  _  ^0 

w0  % 


2.73 


Prom  equation  2.59,  the  unnormalized  frequency  transforma¬ 
tion  will  be 


P 


2.74 


or 


for 


3  =  0 


S  2+00q  2 


SO) 


0 


2.75 


In  the  frequency  domain  the  band-stop  filter  may 
be  obtained  by  taking  reciprocal  of  the  transformation  for 
the  band-pass  filter: 


0 


P 


s  2  +  w 


2 


2.76 


:  nsern  oi'iJtem 


0* 


vU) 


=  q 


V+£b 


\[£  L:'i  qocTs-  >nsd  srict  nx^mob  v;yn  (p9T:l  srIJ  ril 
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From  equation  2.73  and  equation  2.75 


P 


n 


_  s2+V 

“o  "  sw0 


2.77 


or 


S  2  +  CO 

c 


s  +  cOq  2  =  0 


2.78 


where  go  =  u>,  -  con 

c  n  Jc 

The  roots  of  equation  2.78  will  be  poles  of  the  band-pass 
filter  obtained  from  the  poles  of  the  low-pass  filter: 

S  (O  GO  2 

*  {  =  ~4  +  /  (-i)  -  CO  (0  2.79 

S  2  2  i  2 

From  each  of  the  n  poles  of  a  low-pass  filter  there  will 
be  generated  2n  poles,  in  complex  conjugate  pairs,  in 
the  s-plane. 

The  transfer  function  for  a  normalized  band-pass 
filter  can  be  obtained  using  equation  2.59  in  equation 
2.36: 


YNB(s) 


(s2+l-p  s)  (s2+l-p  s)  (s2+l-p  s)  (s2+l-p  s) 

1 _ 2  3  4 


2.80 


« GO  -  ,U)  = 

n  o 

3£  c  Bd  eriJ  *io  seloq  9d  II iw  8Y-S  .  >1:  .-.fj  i  d 

r'lsllxl  33Bq-woI  9rid  1o  s^Ioq  9ud  rro'i'i  -  9jTj  J 

ro) 

ev.s  .ui  u  -  (§-)  \  +  —  =  )  , 

■ 

8S£>,  i  b9sIl6fmon  b  io'i  ioltfonul  'is'lsas'id  srfT 


(3  q-I+s s )  (3  q-I+'s)  (3  q-I+sa)  (a  q~I+se) 
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or 


YNB(s  ^ 


S8+g  S7+g  s  6+g  s  5+g  s  **  +g  s3+g  s2  +  g  s  +  g 

7  6  5  4  3  2  1  Q 


2.81 


From  equation  2.79  the  transfer  function  for  a  unnormal¬ 
ized  band-pass  filter  can  be  found  as 

s  s  s  s 

Yb(s)  =  - - -  •  -  •  -  •  - 

(s2+b  s+c  )  (s2+b  s+c  )  (s2+b  s+c  )  (s2+b  s+c  ) 

1  1  2  2  3  3  4  4 

2.82 

The  coefficients  b.  and  c.  are  obtained  from  the  pole 

J  J 

positions  which  occur  in  complex  conjugate  pairs  to  make 
each  of  the  terms  in  the  bracket  of  equation  2.82: 

(s  — s  . )  (s  — s*  . )  =  (s2+b.s+c.)  2.83 

J  J  J  J 

The  s-plane  representation  of  equation  2.82  must  be 
converted  to  the  digital  or  z-plane  representation  by  a 
bilinear  z-transform: 


s  2+b . s  ,+c  . 
d  j  d  j 


2At  1-z 

( At ) 2  ( 1+z ) 2 


(1+z) 


4  .  (l-z)2  +  2b,1  ~At  _  (l-z) 

(At)2  (1+z)2  (At)2  (1+z)2 


( 1+z ) +c . 

J 


(1+z) 2 (At) 2 
( 1+z ) 2 ( At )  2 


:S8.S  noidsjjps  lo  ^9?foBrxcl  9iid  nl  srraatf  9 rid1  lo  does 


JA.  ,dS, 


!  (s+I) 


(S+I) 
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~  (i-z2) 

h 


1  +  a. (cj At"At } 
J 


*1(_2  _ 

a .  At 

J 


b  .  +c  . 
J  <3 


At  s 
z  ' 


2.84 


where 


TT  +  b  . 
At  j 


n  A  t 


2.85 


Equation  2.84  can  be  written  as 


where 


+b  .  s+c  . 
J  J 


Ii(1"z2) 

J 


2.86 


B.  =  1  -  FI .  •  z  +  F2 . •  z  2  2.87 

J  J  J 


FI.  and  F2 .  are  the  coefficients  of  z  and  z2  in 
J  J 

the  equation  2.84.  Omitting  the  static  gain  factor 
- - - — ,  the  z-transform  of  a  Butterworth  band-pass 

3  .cl  *3  •  3 

12  3  4 

filter  with  8  poles  is  given  by 


(1-z2 ) 4 

W(z)  =  -  2.88 

B  (z)  B  ( z ) • B  ( z )  •  B  (z) 

12  3  4 


In  the  following  example  a  band-pass  filter  with 
8  poles 3  cut-off  frequences  1  and  10  c.p.s.  and  a  digiti¬ 
zing  interval  of  1/118  sec.  will  be  calculated.  First 
the  cut-off  frequencies  are  corrected  using  equation  2.21 
and  from  equation  2.79  the  roots  are  calculated . (Figure  2.13) 
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id) 

f  (c.p 

-47.0526 

±125.8533 

4.03 

-20.1726 

±129.7337 

9.51 

-  6.6029 

±i  3.6280 

0.58 

-  2.0527 

±i  6.0782 

0.97 

Figure  2.13.  Poles  and  zeros  of  a  Butterworth  band-pass 

filter  with  cut-off  frequencies  at  1  and  10 
cycle  per  second. 
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The  coefficients  of  the  polynomials  in  equation  2.81 


gQ  =  2.67661  x  10 1 0 
g  =  1.00427  x  1010 

1 

g  =  2.14872  x  109 

2 

g  =  2.81529  x  108 


g  =  2.16737  x  10 

4 

g  =  6.96027  x  10 

5 

g  =  1.31337  X  10 

6 

g  =  1.51762  x  10 


The  coefficients  for  equation  2.82  are 


b  =  94.10519 

1 

b  =  40.34528 

2 

b  =  13.20584 

3 

b  =  4.10533 


=  2882.33984 
=  3975.04687 
=  56.76093 
=  41.15739 


The.  poles  of  equation  2.88  are  (Figure  2.14) 


z  =  1.452137  ±  10.335533 

1 

z  =  1.031333  ±  10. 562206 

2 

z  =  1.057052  ±  iO. 032561 

3 

z  =  1.016185  ±  iO. 052410 

4 

The  terms  in  the  denominator  of  equation  2.88  are 

B  =  1  -  1.307476  z  +  0.450190  z2 
1 

B  =  1  -  1.494986  z  +  0.724783  z2 


are 
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B  =  1  -  1.890261  z  +  0.894119  z2 
3 

B  =  1  -  1.962924  z  +  0.965829  z2 

The  impulse,  amplitude  and  phase  response  of  this  filter 
are  shown  in  figure  2.15,  2.16  and  2.17. 


Figure  2 . 14 . 


Location  of  poles  and  zeros  of  equation  2.88. 
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Figure  2.15.  Zero  phase  impulse  response  of  Butterworth 

band-pass  filter  with  a  cut-off  at  1  c.p.s. 
and  10.0  c.p.s.,  for  n=4  and  At=l/ll8  sec. 
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Figure  2.16.  Zero  phase  amplitude  response  of  Butterworth 

band-pass  filter  with  a  cut-off  frequencies 
at  1  c.p.s.  and  10.0  c.p.s.  for  n=4  and 
At=l/ll8  sec. 
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Figure  2.17.  Phase  response  of  Butterworth  band-pass 

filter  with  a  cut-off  frequencies  at  1  c.p.s. 
and  10  c.p.s.  for  n=4  and  At=l/ll8  sec. 
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2.5  Zero  phase  filters. 

The  phase  response  of  a  Butterworth  filter  as 
calculated  above  does  not  have  a  linear  phase  spectrum. 
According  to  equation  2.10  it  would  be  desirable  to  have 
a  filter  with  a  linear  or  zero  phase  shift  so  that  differ¬ 
ent  frequency  components  do  not  have  differential  lags. 
There  are  two  methods  which  are  applicable  to  a  finite 
length  of  data  to  have  a  zero  phase  response: 

1.  The  input  data  is  convolved  with  an  impulse 
response  F(z)  to  produce  an  initial  output.  This  out¬ 
put  is  then  reversed  in  time  and  convolved  again  with  the 
same  operator.  The  second  output  is  reversed  in  time  to 
obtain  the  final  output  signal. 

A  zero  phase  shift  Butterworth  filter  was  cal¬ 
culated  and  its  impulse  response  obtained  for  a  four 
second  interval  of  time.  Since  the  impulse  response  is 
symmetric  about  the  time  of  occurrence  of  the  impulse,  we 
can  say  that  the  filter  is  zero-phase.  The  phase  spectrum 
as  calculated  by  taking  a  fast  Fourier  transform  of  the 
impulse  response  was  zero  to  five  significant  figures.  The 
amplitude  response  had  a  sharper  cut-off  because  two  con¬ 
volutions  were  involved  so  that  the  amplitude  spectrum  is 
the  square  of  the  spectrum,  F(z),  of  a  single  Butterworth 
filter.  That  is,  if  we  let  o(t)  be  the  output  from  a 
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convolution  of  the  input  data  i(t),  and  the  impulse  res¬ 
ponse,  f(t),  then  the  first  output  is: 

oi (t)  =  i( t ) *  f (t) 

The  reversed  output  is: 

o  (t)  =  o  (-t)  =  i(-t)  *  f(-t) 

2  1 

Then  data  is  filtered  once  more  to  give: 


o  (t)  =  {  i ( -t )  *  f(-t)  }  *  f(t) 

3 


The  final  output  involves  a  time  reversal  which  yields: 


o  (t)  =  {  i(t)  *  f(t)  }  *  f(-t) 

4 


Taking  the  Fourier  transform  of  the  last  output 


0  (f)  =  {  1(f)  •  F(f)  }  •  F  (f) 

4 


=  1(f)  •  {  F(f )  •  F  (f)  } 


=  1(f)  •  | F( f ) | 


2 


2.89 


el  duqdao  b seievs*!  SflT 


(3-)l  *  {  (; )1  *  (J)i  }  =  (J)  o 
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where  F  (f)  is  the  complex  conjugate  of  the  transfer  function. 
For  filters  in  which  the  shape  of  the  filter  in  the  pass- 
band  is  important,  it  is  necessary  to  design  F(z)  to  have 
an  amplitude  spectrum  that  is  the  square  root  of  the  de¬ 
sired  spectrum. 

To  process  field  records  with  a  zero  phase  filter 
a  sufficient  number  of  zeros  must  be  added  at  the  end  of 
the  data  to  obtain  a  reasonably  small  amplitude  before  the 
reversed  data  is  filtered  again.  Generally,  between  100 
and  300  zeros  were  found  to  be  sufficient  so  that  the  im¬ 
pulse  response  became  negligible.  The  addition  of  too  many 
zeros  may  give  underflow  problems.  Therefore  both  the  im¬ 
pulse  response  and  the  transfer  function  should  be  studied 
before  a  particular  filter  is  used  on  any  problem. 

2.  Shanks  (1967)  has  advocated  the  use  of  a 
second  method  of  obtaining  a  zero  phase  shift  filter. 

First  a  stable  positive-time  and  a  stable  negative-time  fil¬ 
ter  are  calculated.  The  negative-time  filter  is  obtained 
by  replacing  z  with  1/z  in  the  theory  developed  previous¬ 
ly.  The  zero  phase  output  of  this  filter  can  be  computed  by 
filtering  the  input  with  the  forward  filter  and  adding  to 
this  result  the  input  filtered  backwards  with  the  negative¬ 
time  filter.  Time  reversal  of  the  data  is  obtained  by 
replacing  z  in  F(z)  with  1/z  .  Since  z 
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a  time  advance,  using  equation  2.17,  a  forward  recursion 
relation  can  be  expressed  as 


y 


n 


N 

T  a. x  ,  . 
i=0  1  n+1 


M 

j=lbJ'yn+J 


2.90 


n=k ,  k-1 ,  k-2 , . . . . 


The  zero  phase  impulse  response  was  calculated 
for  a  band-pass  Butterworth  operator,  it  has  a  symmetric 
impulse  response  about  the  time  of  the  input  impulse  but 
the  function  is  entirely  different  from  the  one  calculated 
by  method  1.  The  amplitude  and  phase  spectra  may  be  cal¬ 
culated  using  the  fast-Fourier  transform.  It  was  found 
that  it  had  a  zero  phase  spectrum  but  the  amplitude  spectrum 
is  distorted  in  an  undesirable  manner. 

Steps  followed  in  this  technique  can  be  formula¬ 
ted  as  follows : 

The  output  from  the  negative  recursion  filter 

o  (t)  =  i(t )  *  f (t) 
i 

The  output  from  the  forward  or  negative  time  recursion 


filter  is 
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o  (t)  =  i(t)  *  f ( -t ) 
2 


The  final  output  is  a  sum  of  above  two  equations 


o  (t)  =  o  (t)  +  o  (t)  =  i ( t )  #  {  f(t)+f(-t)  } 

3  1  2 

Taking  a  Fourier  transform  of  last  output  we  have 

0  (f)  =  1(f)  •  (F(f)+F*(f)  }  2.71 

3 


If  we  let  the  real  and  imaginary  part  of  F(f) 
be  FR(f)  and  Fj(f)  ,  then 


0  (f)  =  1(f)  {  2  Fr  } 

3  n 

whereas  the  required  filter  should  have  an  amplitude 
response  resembling 


| F(f ) |  =  /  F  2  +  F  2 
n  1 

Because  the  distortion  due  to  the  loss  of  F^ 2  is  unpredic¬ 
table  for  the  Butterworth  formulation  this  second  method 
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2.6  Conclusion. 

The  Butterworth  recursion  filter  is  a  best 

approximation  to  an  ideal  filter  and  in  computation  the  re¬ 
cursive  form  is  quite  efficient  and  fast.  The  program 
listed  in  Appendix  3  is  able  to  filter  a  6  x  6300  pt. 
array  in  less  than  2.27  minutes,  but  this  time  includes  cal¬ 
culation  of  the  filter  coefficients.  One  can  reduce  this 
time  considerably  by  calculating  filter  coefficients  by  a 
separate  program.  It  was  noticed  that  the  calculation  speed 
was  increased  by  the  number  of  records  in  the  IBM  360/67 
system.  For  example  3  records,  each  approximately  51 
seconds  long,  took  about  5  minutes  to  filter  with  the  same 
program.  The  transfer  function  and  impulse  response  for 
one  filter  which  has  proven  to  be  very  useful  is  shown 
in  figures  2.18  and  2.19.  Examples  of  unfiltered  records 
are  shown  in  figure  2.20  while  the  same  data  after  digital 
filtering  is  shown  in  figure  2.21  and  2.22.  It  is  obvious 
that  the  first  and  later  arrivals  on  the  records  are  inter¬ 
preted  much  more  easily. 
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Figure  2.18.  Impulse  response  of  zero  phase  shift 

Butterworth  filter  for  cut-off  frequencies 
at  1  c.p.s.  and  5  c.p.s.  The  sampling 
interval  is  1/118  sec. 
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Figure  2.19.  Amplitude  response  of  zero  phase  shift 

Butterworth  filter  for  cut-off  frequencies 
at  1  c.p.s.  and  5  c.p.s.  The  sampling 
interval  is  1/118  sec. 


-^lon9upe  :  1';  riol  -  *  tl  •  1  nJiovn  -^d^u" 


AMPLITUDE  RESPONSE  OF  BUTTERWORTH  BAND  PASS  FILTER  FOR  1.0-  5.0  HZ 
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Figure  2.20.  Example  of  4  unfiltered  seismic  records  from 

Project  Early  Rise 
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Figure  2.21.  Four  seismic  records  from  Project  Early  Rise 

filtered  by  a  zero  phase  shift  band  pass 
recursive  digital  Butterworth  filter  with  cut¬ 
off  frequencies  at  1  c.p.s.  and  5  c.p.s. 
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Figure  2.22.  Four  seismic  records  from  Project  Early 

Rise  filtered  by  a  zero  phase  shift  band 
pass  recursive  digital  Butterworth  filter 
with  cut-off  frequencies  at  0.5  c.p.s.  and 
3.0  c.p.s. 
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CHAPTER  3 

THE  CRUSTAL  TRANSFER  FUNCTION 

AND 

SPECTRAL  RATIOS 


3.1.  Theory. 

A  solution  of  the  boundary  value  problem  for 
seismic  waves  transmitted  in  a  layered  medium  by  the 
Thomson-Haskell  matrix  formulation  gives  a  rapid  and  exact 
calculation  of  the  corresponding  transfer  function.  Bound¬ 
ary  conditions  for  the  solution  of  the  problem  are  the 
continuity  of  stress  and  displacement  at  each  internal  in¬ 
terface  and  vanishing  of  stresses  at  the  free  surface.  It 
is  assumed  that  for  the  earth's  crust  the  ratio  of  longitu¬ 
dinal  to  shear  wave  velocity  is  approximately  constant  (/3), 
corresponding  to  a  value  of  0.25  for  Poisson’s  ratio.  It 
is  also  assumed  that  the  density  increases  linearly  with 
increasing  seismic  velocities  (Birch,  1964). 

The  theory  for  obtaining  the  transfer  function 

has  been  given  by  Haskell  (1962),  (see  Appendix  1).  The 

transfer  functions  for  the  horizontal  and  vertical  components 

of  ground  motion  at  a  given  frequency  for  a  plane  P-waves 

t  h 

which  is  incident  at  the  bottom  of  the  n  layer  are  given 
by 
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where  c  is  the  apparent  surface  velocity, 


c>a 
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D  =  (J  -  J  )(J  -  J  )  -  (J  -  J  )  ( J  -  J  ) 

11  21  32  42  12  22  31  41 


3.3 


and  Jpp  are  the  elements  in  the  matrix 


- 1 


J  =  E._  *a  ,  .  a  ~  .  a. 


n 


♦  a  .  a.  0 
n-1  n-2 


3.4 


E"1  is  a  4x4  matrix  associated  with  the  boundary 

condition  between  the  half-space  and  the  n-1^ 
layer.  This  matrix  is  independent  of  frequency. 


a  t  is  a  4x4  matrix  associated  with  transmission 
n-1 

across  the  n-l^h  layer  and  the  boundary 
between  the  n-lth  layer  and  n-2n(^ 
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This  matrix  depends  on  the  layer  constants, 
the  angle  of  incidence  and  the  frequency  of  the 
input  plane  wave. 
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Figure  3.1.  Notation  used  for  a  crustal  system  of  n  layers. 

Similarly  the  horizontal  and  vertical  component 
of  ground  motion  for  an  incident  S  wave  can  be  written  as 
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A  fortran  IV  program  written  by  Hanon  to  calculate 
the  transfer  function  has  been  modified  for  operation  on 
the  University  of  Alberta's  IBM  360/67  (Appendix  2). 

3.2.  Properties  of  the  Transfer  Function. 

In  order  to  investigate  the  crustal  structure 
using  the  transfer  function  it  is  necessary  to  know  how  it 
behaves  with  variations  in  the  crustal  parameters.  Hanon 
(1964)  and  Fernandez  (1965)  have  studied  the  behavior  of 
the  crustal  transfer  function.  It  will  be  useful  to  review 
briefly  some  of  Hanon' s  conclusions: 

1.  "Variations  in  the  total  thickness  of  the 
crust  as  small  as  one  kilometer  cause  a  noticeable 
shift  in  the  position  of  the  peaks  and  troughs  of 
the  transmission  coefficients." 
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2.  "For  crustal  models  having  the  same  total 
thickness,  the  main  difference  in  the  frequency  varia¬ 
tion  of  the  transmission  coefficients  is  in  the  charac¬ 
ter  of  the  os cillat ions . " 

3.  "The  variation  in  the  horizontal  transmission 
coefficient  is  more  irregular  than  that  of  vertical 
coefficient  in  the  range  of  angles  of  incidence  studied." 

In  this  study  a  crustal  model  (Table  1.)  obtained 
by  Cumming  and  Kanasewich  (1966)  for  southern  Alberta  is 
used . 


TABLE  I 

-  CRUSTAL  MODEL 

FOR  SOUTHERN 

ALBERTA 

Depth  (km) 

Vp(km/sec ) 

Vg( km/sec ) 

Density  (gm/ 

1.105 

2.31 

1.33 

2 . 04 

2 . 202 

3.06 

1.77 

2.21 

33.4 

6.50 

3.75 

2.73 

43.2 

7.15 

4.13 

3.20 

oo 

o 

oo 

4.66 

3.45 

The  transfer  function  and  spectral  ratios  will  be  calcula¬ 
ted  for  this  model. 

When  the  thickness  of  the  third  layer  was  increas¬ 
ed  by  3  kilometers  this  caused  noticeable  shifts  in  the 
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peaks  of  the  transfer  function  (Figure  3.2).  The  sensiti¬ 
vity  of  the  transfer  function  to  variations  in  the  P-wave 
velocity  in  the  same  layer  was  also  tested  by  increasing 
the  velocity  by  1  km/sec.  This  introduced  a  change  of 
amplitude  in  the  transfer  function  as  well  as  shift  in  the 
peaks  (Figure  3.3).  The  above  test  is  also  performed  for 
a  different  angle  of  incidence  and  as  seen  in  figure  3.^ 
it  is  insensitive  to  variations  in  this  parameter.  The  phase 
response  of  the  transfer  function  of  above  model  is  almost 
linear  as  seen  in  figure  3.5. 

3.3.  Application  of  the  Crustal  transfer  function. 

A  signal  transmitted  through  the  crust  and  re¬ 
corded  at  the  surface  is  modified  in  a  manner  depending  on 
the  density  and  elastic  properties  of  the  media.  The 
spectrum  of  such  a  signal  can  be  expressed  as  a  convolution 
of  a  source  function  with  the  crustal  transfer  function: 

R ( f )  =  P(f)  *  H ( f )  3.7 

where  P(f)  is  the  spectrum  of  the  incident  pulse  and 
H ( f )  is  the  transfer  function  which  depends  upon  the 
crustal  parameters.  Equation  3.7  can  be  written  for  both 
the  vertical  and  horizontal  components: 
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Figure  3.2.  Modulus  of  transfer  functions  and  their  ratios 

for  plane  waves  incident  on  a  southern  Alberta 
crustal  section  at  an  angle  of  incidence  of 
80°.  The  depth  has  been  increased  in  the  third 
layer  by  3  kilometers.  The  crosses  denote  the 


thicker  model 


VERTICAL  TRANSFER  FUNCTION 


64-A 


FREQ 


65 


Figure  3.3.  Modulus  of  the  transfer  function,  their  ratios 

and  the  phase  of  the  ratios  for  plane  waves 
incident  on  a  southern  Alberta  model  at  an 
angle  of  incidence  of  80°.  The  velocity  has 
been  increased  in  the  third  layer  from  6.50 
to  7.50  km/sec.  (see  table  1).  The  crosses 
denote  the  higher  velocity  model 
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Figure  3.^.  Modulus  of  transfer  function  and  their  ratios 

for  plane  waves  incident  on  a  southern  Alberta 
crustal  model  at  an  angle  of  incidence  of  60°. 
The  depth  has  been  increased  in  the  third  layer 
by  3  kilometers.  The  crosses  denote  the  thicker 


model 
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Figure  3.5.  Phase  curves  for  transfer  functions  of  a 

southern  Alberta  crustal  model.  This 
corresponds  to  the  modulus  of  transfer  func¬ 
tion  shown  in  figure  3.2. 
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Rz(f)  =  P(f)  *  Hz(f) 


3.8 


RH(f)  =  P(f)  *  HR(f) 


3.9 


The  crustal  transfer  function  can  be  calculated  and  if 
one  knows  the  source  function,  P(f),  it  is  possible  to 
obtain  the  recorded  signal  using  equation  3.7.  Since  the 
impulse,  P(f),  is  not  known,  it  is  usual  to  take  the  ratio 
of  the  horizontal  and  vertical  spectrum: 


Rz(f)  Hz(f) 


3.10 


RH(f)  HR(f) 


By  taking  the  ratio,  the  unknown  function  P(f)  is  elimina¬ 
ted.  The  observed  spectral  ratio  can  be  matched  to  the 
theoretically  calculated  ratio  of  transfer  functions  for 
any  crustal  model. 


Phinney  (1964)  was  the  first  to  use  the  ratio  of 


the  transfer  functions  to  obtain  crustal  structure  using 
long  period  data.  Fernandez  (1965)  calculated  master  curves 
to  determine  the  crustal  structure  from  the  spectrum  of 
long  period  P  waves  in  the  United  States.  In  their  studies 
a  time  window  applied  to  observed  spectra  was  neglected  be¬ 
cause  long  period  data  was  used.  For  short  period  data  the 
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window  must  be  short  so  that  additional  phases  do  not  arrive 
and  its  effect  cannot  be  neglected  (Leblanc,  1966).  A 
recorded  seismogram  seen  through  a  window  W(f)  can  be 
expressed  as : 


RT(f)  =  R(f)  »  W(f)  3.11 

or  in  the  time  domain 

rT(t)  =  r ( t )  •  W(t)  =  R(f)  *  W(f)  3.12 

using  equation  3.7  in  equation  3.11: 

Rip(  f )  =  [P  ( f )  •  H(  f )  ]  *  W(f)  3.13 

A  difficulty  arises  in  equation  3.13  because  the  convolu¬ 
tion  must  be  performed  on  two  unknown  functions.  However 
Leblanc  (1966)  showed  that  the  convolution  can  be  done  on 
one  unknown  only;  e.g.  on  H(f): 

Rrp(f)  =  [P(f)-H(f)]  *  W(f)  =  P(f)  •  [H( f ) *W( f ) ] 

T  3.14 

This  convolution  is  valid  under  two  conditions  (Leblanc, 

1966 )  : 
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1.  The  transfer  function  in  the  time  domain 
can  be  approximated  by  a  sequence  of  weighted  delayed 
delta  functions. 

2.  The  source  function  is  assumed  to  be  of 
finite  duration  and  short  compared  to  time  duration 
of  the  record  T  ,  to  prevent  truncation  of  any  sub¬ 
pulses  by  the  window.  Therefore  equation  3.10  can  be 
written  as 


RZT(f)  _  HZT(f) 
RHT(f)  HRT(f) 


3.15 


A  convolution  of  the  transfer  function  with  a 
window  can  be  performed  by  a  method  explained  by  Leblanc, 
1966,  but  there  are  difficulties  in  the  case  of  a  rectangu¬ 
lar  window  as  indicated  because  of  its  slow  convergence. 

In  this  study  the  truncation  of  the  transfer  function  has 
been  done  by  a  double  Fourier  transform.  A  simple  rectangu¬ 
lar  window  W(t)  is  used  where 

W(  t )  =  1  o<t<T 


o 


elsewhere 
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and  T  is  the  time  duration  of  the  signal.  The  steps  used 
in  truncation  are: 

1.  The  transfer  function  was  calculated  using 
the  program  listed  in  Appendix  2.  The  time  synthesis 
of  transfer  function  was  obtained  by  a  Fourier  trans¬ 
formation  using  a  Fast  Fourier  algorithm  (Cooley  and 
Tukey ,  1965 )  : 


-2iriWT 

N 


h(t)  =  i  l  H(f)  e 
T=o 


3.17 


In  the  expression  above  N  is  the  total  number  of  points 
calculated  at  intervals  of  Af  c.p.s.  and  W  is  the 
digital  frequency  index. 

2.  The  desired  length  of  h(t)  is  generally 
chosen  by  a  study  of  time  synthesis  of  the  transfer 
function.  This  time  length  must  be  equal  or  almost 
equal  to  the  time  duration  of  the  seismic  records 
under  analysis  but  at  the  same  time  such  that  h(t) 
has  a  value  as  near  to  zero  as  possible  for  t  =  T. 
Taking  the  Fourier  transform  of  equation  3.17 


N-l 

H ( f )  =  l  h(t)  e 
T=o 


-27TiWT 

N 


3.18 
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Figure  3.6  shows  the  steps  described  above.  Here 
4.3^  seconds  window  was  used  since  the  major  pulses  have 
arrived  at  this  time.  A  perfect  time  synthesis  of  the  trans¬ 
fer  function  would  be  obtained  using  the  full  range  of  inte¬ 
gration  i.e.j  from  -°°  to  +°°  .  For  a  narrower  range  of 
frequency  one  must  assure  that  the  transfer  function  is  not 
affected  by  the  rejected  frequencies.  For  this  reason  the 
stability  of  the  truncated  transfer  function  was  tested  with 
calculation  using  higher  frequencies  in  the  time  synthesis. 

It  was  noted  that  when  the  synthetic  seismogram  was  trunca¬ 
ted  over  a  constant  time  interval  (4.34  seconds)  that  the 
truncated  transfer  function  in  the  frequency  range  0-10 
c.p.s.  was  similar  to  the  truncated  transfer  function  with 
the  frequency  range  0  -  15  c.p.s. 

3.4.  Analysis  of  data. 

Power  spectra  of  field  records  were  computed  by 
first  calculating  the  autocorrelation  function,  this  func¬ 
tion  was  modified  by  a  Daniell  window  whose  length  is  4.34 
seconds  and  the  power  spectrum  was  obtained  by  a  fast 
Fourier  transform  of  the  modified  autocorrelation  function. 
The  field  records  were  filtered  with  a  Butterworth  band¬ 
pass  filter  modified  for  zero  phase  shift,  whose  cut-off 
frequencies  are  1  and  5  c.p.s.  This  filter  reduced  the 
noise  background  on  the  power  spectrum  (Figure  3.7).  The 
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Figure  3.6.  Truncation  of  transfer  function  for  Alberta 

model  using  frequency  interval  of  Af=0.05 
and  time  interval  of  At=l/50 
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main  peak  is  centered  at  about  2  c.p.s.  ,  small  peaks  at 
higher  frequencies  on  some  records  appear  to  represent 
mainly  noise  since  the  peaks  are  not  consistent  from  trace 
to  trace  and  record  to  record. 

For  Project  Early  Rise  the  vertical,  horizontal 
and  ratio  spectrums  were  calculated  at  9  locations  in  Alberta 
and  Saskatchewan.  Theoretical  spectral  ratios  were  also 
calculated  for  a  southern  Alberta  model  (Table  1)  to  compare 
with  the  field  records.  At  location  WC  25,  which  is  the 
beginning  of  the  line  recorded  by  the  University  of  Alberta, 
the  theoretical  spectral  ratio  obtained  for  this  model  agrees 
well  with  the  spectral  ratios  from  the  field  records  (Figure 
3.7).  In  this  comparison  one  must  make  use  of  only  the 
part  around  the  mid-frequency  rather  than  the  higher  or  lower 
frequencies  which  are  distorted  by  various  instrumental 
effects.  The  last  recording  station  in  Saskatchewan,  WC  34, 
did  not  agree  with  these  theoretical  spectral  ratios  and 
an  attempt  has  been  made  to  find  a  better  model  by  changing 
the  depths  of  the  layers.  After  11  trials,  it  was  found 
that  the  spectral  ratio  for  a  model  with  the  third  layer 
thickness  increased  to  about  36  kilometers,  fits  the  spectral 
ratios  of  the  field  record  better  (Figure  3.8).  In  Alberta 
the  spectral  ratios  for  the  field  record  at  WC  39  was  found 
to  be  similar  to  the  theoretical  spectral  ratios  when  the 
third  layer  thickness  was  increased  to  about  4l  kilometers 


■ 

.  3  -:cD9'i  bLsl.  :  3rict  i  /  w 

j  v,  n  )  Jo  e*u  a  10  raa.i-j-  ,  ".do  •- .if  nT 

'  9W0.  '  '  -ceris  tr.  :d  -Id  •iorldfii  ye  j  9  ":V-5.;  -  :  t  n  •  >  ".b  --£<1 

iBtfne:  iirfeni  s  /c  .  isv  v,d  b.i.  co:  ci  3  nolfw  a  'or  -  T 

9rid  neriw  soidB-i  ifi'idosqs  LBOldeia^dS  srid  c  d  rcBli  rtie  sd  od 
3,i9demolJi.3l  :  duodB  od  b93B9iiont  8«v/  as  rDioirid  ri9\  I 


75 


(Figure  3.9). 

The  spectral  ratios  for  the  field  records  were  not 
very  stable  from  one  location  to  the  other.  The  change  in 
the  spectral  behavior  can  be  explained  by  three  effects: 

(1)  variation  in  the  source  function,  (2)  upper  mantle 
effects  due  to  triplication  in  the  travel  time  curve  and 
(3)  crustal  effects  such  as  non  parallel  layering  or  fault¬ 
ing.  It  is  difficult  to  predict  the  pulse  shape,  however, 
the  pulse  can  be  approximated  by  various  functions  such  as 
the  Sin  X/X  which  was  used.  Since  the  charge  was  detona¬ 
ted  in  the  water  at  a  location  which  varied  slightly  from 
shot  to  shot  it  has  been  impossible  to  represent  these  varia¬ 
tions  in  the  pulse  function.  However  the  effect  of  the  pulse 
can  be  eliminated  to  great  extent  by  taking  spectral  ratios. 
One  must  also  remember  that  the  Lake  Superior  crustal  struc¬ 
ture  where  the  source  was  generated  can  effect  the  spectrum 
of  field  records.  The  most  likely  explanation  for  the  rapid 
changes  in  the  character  of  the  spectrum  is  the  effect  of 
the  upper  mantle.  Body  waves  traveling  from  the  shot  point 
to  the  recording  station  are  modified  by  mantle  structure 
and  in  particular  reflected  waves  are  generated  at  second 
order  discontinuities.  Interference  of  one  wavefront  with 
another  can  produce  minima  and  maxima.  The  observed  spectrum 
was  found  to  agree  with  the  theoretical  one  obtained  for  a 
reasonable  Earth  model  and  seems  that  the  technique  can  be 


•  •  os  ■  ri  t  i  t  ~ .  i  x €  x 1  i i  T"  :  d  ..  '■  -  ■  :  j 


8G  a  s  rjfJDrwJ  suotiisv  yd  mixo«iqqfi  i  n;o.  '  uq  ■ 
-  6  asv.  .  :£,ii'  erfj  DnJ. 


r  1  Ot'iiBv  rfo Jt  v;  no.-: isool  8  ds  d  nl  be! 

.eoxJ'  .  J- o  c2  ■  nxxfid  yd  nsdx  d£  r  o  l-  o 


76 


successful  if  used  with  discrimination. 

3.5.  Conclusion. 

The  transfer  function  calculated  by  the  Thomson- 
Haskell  matrix  method  has  been  found  to  be  very  sensitive 
to  the  variation  of  crustal  parameters.  Spectral  ratios 
for  field  records  were  compared  with  the  theoretical  ones 
with  the  thicknesses  of  the  layers  being  varied.  It  appears 
likely  that  the  crustal  thickness  increases  in  the  west¬ 
ward  direction  along  the  line  of  recording  stations.  How¬ 
ever  for  a  limited  number  of  recordings,  we  could  not  find 
a  plausible  crustal  model.  One  can  get  a  better  fit  of  the 
theoretical  spectral  ratios  to  the  observed  spectral  ratios 
by  changing  parameters  other  than  thickness  but  the  amount 
of  model  computation  can  be  very  large.  In  general,  it 
appears  that  the  spectral  ratios  of  short  period  body  waves 
is  very  sensitive  diagnostic  for  obtaining  crustal  structure 
if  one  can  obtain  a  reasonable  level  of  signal  to  noise  from 
sources  at  distances  where  the  upper  mantle  complications 
do  not  affect  the  crustal  reverberations.  It  is  also 
difficult,  at  this  stage  to  determine  if  one  has  a  unique 
model  or  only  one  of  a  family  of  possible  ones.  Recording 
of  spectral  ratios  over  a  wider  band  of  frequencies  could 
resolve  this  ambiguity. 
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Figure  3*7.  Vertical  and  horizontal  spectrum  of  a  field 

record  at  location  WC  25  and  comparison  of 
the  spectral  ratios  with  a  theoretical  model 
which  is  represented  by  a  dashed  line 
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Figure  3.8.  Vertical  and  horizontal  spectrum  of  a  field 

record  at  location  WC  3^  and  a  comparison  of 
spectral  ratios  with  the  theoretical  model 
which  is  represented  by  a  dashed  line 
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Figure  3.9.  Vertical  and  horizontal  spectrum  of  a  field 

record  at  location  WC  39  and  a  comparison  of 
the  spectral  ratios  with  the  theoretical  model 
which  is  represented  by  a  dashed  line 
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CHAPTER  4 

SYNTHETIC  SEISMOGRAMS 


4 . 1  Theory . 

A  time  history  of  the  surface  motion  of  a  layered 
system  responding  to  a  pulse  is  obtained  by  convolving  the 
transfer  function  with  a  function  simulating  the  pulse. 

Hanon  (1964)  first  obtained  a  time  synthesis  for  a  number  of 
plane  layers  using  the  Thomson-Haskell  matrix.  He  studied 
the  effects  of  the  pulse  frequency,  angle  of  incidence, 
and  layering  of  the  system  on  the  seismograms.  McCamy  (1967) 
made  similar  theoretical  calculations  but  when  he  compared 
the  results  with  actual  field  records  he  did  not  find  many 
similarities.  In  this  chapter  synthetic  seismograms  will 
be  calculated  for  a  southern  Alberta  model  of  the  crust  and 
compared  with  field  records  from  Project  Early  Rise. 

4.2  Generation  of  synthetic  seismograms. 

In  this  study  a  simple  pulse  function  (Figure  4.1) 
is  used  to  calculate  seismograms : 


P( f )  =  1  f  <  10  c.p.s. 


4.1 


0 


f  >  10  c.p.s. 
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Figure  4.1. 


The  source  function. 
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The  corresponding  time  function  is 


p(t ) 


_1 

2tt 


00 


p(f) 


ei2irft  _  1_  Sin  20-n-t 
tt  t 


4.2 


A  time  synthesis  of  ground  motion  is  obtained  by  taking  the 
inverse  Fourier  transform  of  the  product  of  the  transfer 
function  of  the  layered  crustal  system,  H(f)  obtained  from 
the  Thomson-Haskell  matrix,  and  the  pulse,  P(f): 

00 

r  ( t )  =  ~  /  H  ( f )  •  P(f)  el2Trft  df  4 . 3 


An  assumption  has  been  made  that  we  have  plane  waves  and 
plane  elastic  layers  in  welded  contact  which  do  not 
attenuate  the  energy  in  any  inelastic  process.  The  seis¬ 
mograms  were  calculated  according  to  the  following  steps 
(also  see  figure  4.2): 

1.  The  modulus  of  the  transfer  function  ,  j  H ( f ) |  ,  and 
the  phase  were  calculated  using  the  program  listed 

in  Appendix  2 . 

2.  The  real  and  imaginary  part  of  the  transfer  func¬ 
tion  were  calculated: 


■ 


9 vi  w  sfiBiii  n99d  e.  ■  1  ;3jb 

jb  i  O  ;  '  .jo  ;ie;:iv>  2  c  :->M  -  -  Q  . 

.  saaoo'iq  r  :_:f  f  L^nx  v-.B  ni  \, 


83 


Figure  4.2. 


Calculation  of  theoretical  seismogram  using 
fast  Fourier  algorithm. 
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Re  H  =  HR  =  |H(f) |  Cos  0 
Im  H  =  Hj  =  | H ( f ) |  Sin  0 

3.  The  real  and  imaginary  parts  of  the  transfer  func¬ 
tion,  H,  were  multiplied  by  the  pulse  function. 

4.  For  fast  Fourier  transformation,  the  real  and 
imaginary  part  of  the  transfer  function  were  made  to 
have  even  and  odd  symmetry  respectively  about  the 
mid-point  or  zero  frequency . 

5.  The  synthetic  seismogram  was  obtained  by  a  fast 
Fourier  algorithm  (Cooley  and  Tukey ,  1965): 

N-l  27riWT 

h ( T)  =  i  l  (Hr(W)  +  i  H-j-CW))  e  N  4.4 

W=o 

W  =  o,l,2  • • • •  N-l 
T  =  o,l,2  • • • •  M-l 

where  M  <  — |  3  W  is  a  digital  frequency  index,  and 
T  is  a  digital  time  index.  N  is  chosen  according 


to  the  relation 
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Figure  4 . 


Theoretical  seismogram  calculated  for 
Hanon’s  (1964)  model 
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N 


1 


At  •  Af 


where  At  and  Af  are  sampling  intervals  in  time 
and  frequency. 

To  check  the  program  using  a  fast  Fourier  trans¬ 
form  a  seismogram  was  calculated  -  using  Hanon ’ s  (1964) 
crustal  model.  A  comparison  (Figure  4.3)  shows  that  we 
can  duplicate  these  results  in  a  satisfactory  manner.  As 
a  further  check  the  same  seismogram  was  calculated  using 
another  Fourier  transform  (Howe,  1958)  with  similar  results. 
Another  test  is  simply  carried  out  for  a  simple  layer  model 
over  an  infinite  half  space  using  different  angle  of  inci¬ 
dence.  The  converted  phases  (Figure  4.4)  and  multiple 
arrivals  may  be  predicted  by  ray  theory  and  these  agree  well 
in  amplitude,  phase  and  time  with  arrivals  on  the  synthetic 
seismogram.  For  this  case  the  difference  in  travel  time 
between  P  arrival  and  S  conversions  is 


At 


c°s  9lp 
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h 
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4.5 
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where  <j>2p  is  the  angle  of  Incidence  below  the  Interface 
for  P  waves. 

4>ip  and  <f)^s  are  the  angles  In  the  layer  for 
P  and  converted  S  waves. 

and  oijj  and  8^  are  P  and  S  wave  velocities  in  the 

o  th  n 

1  layer. 

Prom  Snell’s  law  one  can  also  write: 

a  8  a 

- ! _  =  _ I _  =  _ 2  n  6 

Sin  4>lp  Sin  <|>ls  Sin  <J>2p 

Substituting  equation  4.6  in  equation  4.5  one  can  obtain 

h  a 

At  =  cT  ^ ~  C°S  ^ls  “  Cos  ^lp^  4,7 

i  i  p 

Using  the  assumption,  which  was  made  in  calculations  of 
the  transfer  function,  that  Poisson's  ratio  is  equal  to 

0.25  for  the  crust,  equation  4.7  will  have  the  form  of 

h  _ 

At  =  —  ( Z- 3  C o s  6 1  -  C o s  4> i  )  4.8 

Oc  JL  -L]0 
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Figure  4.4.  Conversion  in  a  single  layer  model. 


After  assuring  that  the  Fortran  IV  program 
designed  to  calculate  synthesized  ground  motion  worked 
properly,  a  theoretical  seismogram  was  calculated  (Figure 
4.5)  for  a  southern  Alberta  crustal  model.  The  seismogram 
calculated  Is  very  complex  because  of  arrivals  of  different 
converted  phases  from  four  Interfaces,  however.  It  is  still 
possible  to  recognize  and  identify  the  late  arriving  conver¬ 
ted  phases o  In  calculating  this  seismogram  the  angle  of 
Incidence  is  taken  to  be  80°,  this  Insures  that  the  differ¬ 
ent  phases  are  separated  sufficiently  in  time  and  is  a 
realistic  approximation  to  field  data  at  distances  of  500 
to  2000  km  from  the  shot.  A  frequency  interval  of 
Af  =  0.02  c.p.s.and  a  time  interval  of  1/50  second  was 
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Figure  *1.5.  Theoretical  ground  motions  for  a  southern 

Alberta  model  using  a  frequency  interval 
Af =0 .02  c  .  p  „  s . 
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Figure  4 0 6 „  Theoretical  ground  motions  for  a  southern 

Alberta  model  using  a  frequency  interval 
A f =0.005  c.p.s. 
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used.  The  same  seismogram  was  calculated  using  a  frequency 
interval  of  Af  =  0,005  c.p.s,  and  a  time  interval  of  1/40 
second  (Figure  4,6),  to  check  the  stability  of  the  seismo¬ 
gram.  A  frequency  interval  of  Af  =  0,02  c.p.s,  was  found 
to  sufficient  to  obtain  a  seismogram  without  any  aliasing 
problems . 

4.3  Field  records  and  theoretical  seismograms. 

Theoretical  seismograms  obtained  above  for  a 
southern  Alberta  crustal  model  (Table  1)  are  compared  with 
field  records.  In  this  comparison,  filtered  field  records 
were  used  since  filtering  enhanced  the  signals  and  removed 
much  of  the  locally  generated  high  frequency  noise  (Figure 
2.20,  2,21,  2.22),  To  make  the  comparison  of  theoretical 
and  observed  data  easier,  the  theoretically  determined 
seismograms  were  also  filtered  with  a  Butterworth  band  pass 
filter  (Figure  4.7), 

The  vertical  and  horizontal  components  of  the 
field  seismograms  show  similarities  with  the  one  obtained  by 
the  Thomson-Haskell  matrix  formula.  Many  field  seismograms 
are  complicated  because  of  arrivals  from  the  mantle.  In 
figure  4.8,  a  theoretical  seismogram  compared  with  the  field 
seismogram  for  two  pulses  of  energy  travelling  by  different 
paths  in  the  upper  mantle.  The  comparison  for  the  two 
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arrivals  show  substantial  agreement  between  the  field 
data  and  the  theoretical  seismograms.  The  comparison  is  not 
very  good  at  certain  locations  (Figure  4.9)  which  means  that 
models  must  be  modified  using  the  results  obtained  in  Chap¬ 
ter  3 . 
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Figure  4.7.  Theoretical  seismograms  for  a  southern 

Alberta  model  and  its  filtered  version  using 
a  Butterworth  band-pass  filter  whose  cut-off 
frequencies  are  1  c.p.s.  and  5  c.p.s. 
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Figure  4 . 
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Figure  4.9.  Additional  comparison  of  vertical  and 

horizontal  components  of  a  field  record  with 
the  theoretical  seismograms  for  a  southern 


Alberta  model 
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APPENDIX  1 


THOMSON-HASKELL  MATRIX 


FORMULATION 


A. 1.1.  Introduction. 


This  appendix  is  the  summary  of  following 


papers:  Haskell  (1953) >  Haskell  (I960),  Haskell  (1962) 
Thomson  (1950).  The  following  notation  is  used  in  this 
section : 

e  =  i!i  +  :  Cubical  dilation 
3x  3y  3  z 


3u. 

— i) 


Rotation  vector  (Anti-symmetric 


vector) 


a  =  /  ^  +  2  -  :  Velocity  of  P  waves. 


X,  y,  p  are  Ldme  constants  and  density. 


$  =  /  —  :  Velocity  of  S  wave. 

P 


k  =  —  =  wave  number 
a  a 


wavelength 


wave  number. 


i  r  ■  ?jr,  r>' 


.noL;  t:h  .'i:  i  I 


‘t  '■  9f,  Tfy ?  :3r  u  -to  tosv  rcoJ  3^0H 


A. 2 

c  =  Apparent  phase  velocity  measured  along  the  surface. 

tu  _  2tt  2tt 

c  c  -  : - 7 — : - - - —  =  T  :  Apparent  wave- 

°  Horizontal  wavelength  A^  ^ 

number  as  measured  along  the  surface. 


'  f 

z 

A. 1.1.  Wavelength  and  wavefront  of  a  wave  propaga¬ 
ting  in  z  direction. 

A 

=  ~  (Figure  Al.l). 

AH 

3u  .  3u . 

=  A0  6.  .  +  y  +  Tw~)  =  stress  tensor. 

1J  dXj 

1  i  =  j  normal  stress 

=  { 


Figure 


Cos  e 


0 


i  5*  j 


shear  stress 


[  J  l  '  ■  -  ■  .  ■  (  :  . 


\ *uj  3  e  •  ;  no^B  b  *x  .  es*  .  *.  n 


,ma  ,ng  =  Dir>ecti°n  cosines  of  the  normal  to  the 

advancing  dilatation  wave  front  with  respect  to 
x  ,y ,  z  axes  . 

=  Direction  cosines  of  the  normal  to  the  rotational 
wave  front  with  respect  to  x,y,z  axes. 


0Q  =  Angle  of  incidence. 


4*^  » 


Amplitude  of  incident,  reflected,  transmitted 
P  waves . 


ft. ,ft  ,ft.  =  Amplitude  of  incident,  reflected,  transmitted 
1  r  u 

S  waves . 


in  an  n 


layer  model. 
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A. 1.2.  Wave  Equations  and  plane  wave  solution. 


A.  4 


Differentiating  cubical  dilatational  with  respect 


to  x : 


30  _  ,3^u 
3x  '  3x2 


11  _  32u 

9  2  v 

9  2  oo 

A. 1.1 

ax  9x2 

9x9y 

9x9  z 

ubtracting  the  term 

c±(32“ 

+ 

i^)]: 

ay2 

9z2 

iii  +  »!«)  + 

9  ^  9v 

9u^ 

9( 

9  u  9  oo  ^ 

A. 1.2 

3y 2  3z2 

9y  9x 

ay 

9  z 

9  z  9x 

Since 


rs  _  9_v  _  9u  d  ^  _  9u  9 (jo 

z  9x  ~  9y  y  9z  “  9x 


And  rearranging  equation  A. 1.2  we  get 


v2u  ■  H  -  2 


957  957 

_ £  +  2  — i 

9y  9  z 


A. 1.3 


Similarly  we  can  write 


n  2  O  U  0 

V^V  =  - —  -  2 


90 

ay 


957 


957 


9  z 


+  2 


9x 


A, 1.4 


V2w  =  ~  -  2 


90 
9  z 


957  957 

£  +  2  x 


9x 


ay 


a. 1.5 


A  plane  dilatational  and  rotational  wave  will  have  the 


form : 


,  •  w  :£i  p^i  r- \ 


w  0gx:j . 


i . 

. 

...  •.  ••  • ..  •  • 

S  V.  1 

.  -j  a..  •-/  ■■  i  j  /  i  *  \  ■  •  j  •»  /«  *  * .  ■  ' 


•  r 


J':i  ■  j:'v 

;  fe>  V  -“Jr  “  U  V 

, 


;  •  ; 


U.I.A 


e.x.A 


"  i :  I  1  .  ......  f.i  ' 


0  =  e 


i{cot-k  (£  x  +  m  y  +  n  z)} 
a  a  cr  a 


A. 


i{a)t-kg(  £^x  +  m^y  +  ngz)} 

—  e 


A. 


Displacement  can  be  written  as 


i{oot-ka(£  x  +  m  y  +  n  z)} 


u0  e 


a  a1 


a 


A. 


Substituting  into  A. 1.3  the  equation  of  wave  will  be 


V2u„  =  - 


k2  ( £ 2  +  m2  +  n2)ufi  =  -  k2  u0  A 
a  a  a  a  0  a  0 


The  solution  of  equation  A. 1.3  can  be  written  as 


u  =  - 


30  _2 


3ft 


z  2 


3ft 


k 


2  3x 


a 


kg2  3y  kB  32 


A. 


Similarly 


v  =  - 


1  30 


0  3ft 
2  x 


3ft 


k  3y  ,2 
a  k6 


3  z 


k 


3x 


3 


A. 


oo  =  - 


1  30 


0  3  ft 
2  _ y 


k2  3Z  ^2  3x 


a 


3 


k 


0  3ft 

_2  _ x 

3y 


3 


A. 


A. 5 

1.6 

1.7 

1.8 

1.9 

1.10 

1.11 

1.12 


It  is  assumed  that  a  plane  wave  is  propagating  in  the 
positive  x  and  z  directions  with  the  amplitude  diminish- 


{  ’>  ;t  +  \  ,'  i  4  x  '  Jw}i 

9  -  0 


Us;  ft  +  V  >  .  *)  1 


r—  -  =  W 

ing  exponentially  in  the  positive  direction.  The  plane 
wave  solution  is  then 


A .  6 


i(a)t-k  a  x)  in  k  z  in  k  z 

0  =  e  a  a  {0 ,  e  a  a  +  0  e  a  a  } 

t  r 


i(a)t-k  £rx)  -in  k  z  inftkRz 

fi  =  e  63  {fi,  e  a  a  +  fi  e  3  3  } 

y  or 


A. 1.13 

A.  1.14 


For  a  two  dimensional  plane  wave  the  cubical  dilatational 
only  components  of  the  rotation  vector  are 


6 


8 u  ,  8  to 

8x  3z 


A. 1. 15 


fi 


y 


1/2  ( 


8u 
8  z 


8co\ 

8x' 


A. 1.16 


By  substituting  in  equations  A. 1.10,  A. 1.11,  A. 1.12  the 
partical  velocities  are  found  to  be 


u  = 


■co£  0 
a 

k 

a 


2tun0  i(tot-k0£Qx) 


3 


3  3' 


k 


3 


{V 


-inekgz 


inftk  z 

fie  6  3  } 
r 


A. 1. 17 


con  i(tot-k  i  x)  -in  k  z  in  k  z  2£R<o 

oj  =  -  —  e  (0te  "  0re  >  +  fiy 


q 


'3 


A. 1.18 


(  X  o  A  ,,  >  -  cfti)  )  1 


Dropping  the  term  e 


A. 7 


A.1.18  and  letting  x  =  0 


icut 


in  equations  A. 1.17  and 


P 


n 


a 


z 


} 


A. 1.19 


Q  =  rigkgZ 


the  equations  become 


0 

u 


oo£  iaj£ 

-j—  Cos  P  {6t+6r}  +  IT2,  Sin  P  {0t“9r}  - 


and 


2ajn 


3 


,  Cos  Q  {fl.-ft  }  + 
Kg  x  r 


i  2con  ( 

k7 


Sin  Q  }  A. 1.20 


—cun  icon 

00  =  — ~  Cos  P  {e^-e„}  +  — 7; — —  Sin  P  { 0 , +0  }  + 


a 


k 


a 


2go£r  i2co£ft 

+  -£-2.  Cos  Q  {^t+fir} - j^-2-  Sin  Q 

3  3 


A. 1. 21 


The  P  wave  velocity  can  be  written  as 


a 


A. 1. 22 


:•  •  ■ 


,  .  '  !. 


I  i;  Si  n  '  '  ’ 


. 


;  •'  <  > 


.  *  u 


I  .  .t  •  ;  . 
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i  .  ;  ..  A 


*  j 


:1  hmt  9is  .-.UxogI-v  r^xr 


3  ,,  / 


v  ■  ,,  ■  .' 


.1  \ 


i  Y  .  '!  V  ;•}  ■ 


ir~ 

'  * *  A  •’  .  .. 

•  ’  \k-  .  \  V.*  •  *  " 


(ft-  +n)  9  800  ~r- 


:  H  ff. 
.0  C- 


ffW 


...  -t 


i"  0  .  ~.v  ?  •  \  <• 

V  -  K  -A  /*•*•  Ui  /  ... 

/*!  A  •?»!•  A.-' 

';jF  j  :::!  . r 

jo 


and  the  horizontal  phase  velocity  as 


A.  8 


c 


A. 1.23 


c 


Figure  A. 1.3 


k 


a 


A. 1.24 


From  figure  A. 1.3  the  following  relations  are  evident: 


£  =  Cos  e  =  77 

a  c 


A. 1. 25 


n  =  Cos  E 
a 


/c2-a2 

c 


A. 1. 26 


r  =  tg  e 
a 


/c2-a2 


a' 


c  >  a 


a 


A. 1.27 


Similarly  S  wave  velocity  can  be  written  as: 


A. 9 


6 


U) 


k 

—  and  c 


0) 

k 


A. 1.28 

A. 1.29 


From  figure  A. 1.4  the  following  relations  also  hold: 


l 


6 


Cos  f  =  — 
c 


A.l. 30 


n 


e 


Cos  F 


/c2-e2 


c 


A.l. 31 


r 


6 


tg  f 


/c  2-g  2 

6 


c  >  6  A. 1.32 


Figure  A. 1.4. 


A.  10 


In  equations  A. 1.20  and  A. 1.21  using  the  relations: 


to 


,  H  Cos (n  k  z)  =  c  —  Cos(kr  z) 
k  a  a  a  2  a 

a  c 


A. 1.33 


and 


2a)nR 

""kT~  Cos(n6  V)= 

P 


r^  Cos(kr^z) 


A. 1.3^ 


It  can  be  written  as 

£  =  -  —  Cos  (kr  z)  {0, +0  }  +  —  Sin  ( kr  z ) { 6  -0  } 

c  2  atr  2  otur 


-  yr  ^  Cos  (kr^z)  +  iyr^  SinCkr^z )  {^t-0>r} 

A. 1.35 

•  2  *2 

—  =  -  —  r  Cos(kr  z){0,-0  }  +  r  Sin(kr  z){0,+0  }  + 

c  2  ot  'a  tr  2  0 1  a  tr 


+  y  Cos  (kr^z)  {$7^+^^}  -  iy  SinCkr^z) A.l.36 

23 2 


where 


Y 


JO 


>  • 


A.  11 


A. 1.3.  Equations  of  stress. 


Normal  stress  can  be  written  as 


1  9  2  0  +  _2  9  2 


9  2  ft 


x 


Pzz  =  Xe  +  2y{-  ^2  3z2  ■  kg  3x3z  -  kg2  3y3z 


}  A. 1.37 


and  shear  stress: 


P  =  2y{ - —  + 

xz 


1_  9^z _ 1_  9  2  _ 1_  329 


2k  2  9x9z  kQ2  9y9z  k  2  9z2  k  2  9z3x 
a  3  3  a 


,  9  2ft 

1 


V  3x1 


.  32n 

+  - —  } 

k  2  9y  9x 


A. 1. 38 


Using  the  equations  A.1.13  and  A.1.14  and  arranging  terms, 
the  equation  A. 1.37  can  be  written  as 


¥  =  -  pa2(y-l)  Cos  P  {8t  +  0r)  +  ipa2(y-l)  Sin  P  {0^.-0^} 


zz 


t  r 


-  pc2y2r^  Cos  Q  {ft^-ft^}  +  ipc2y2rft  Sin  Q  {ft^+ft^} 


't  r 


3 


Jt  r ' 


A. 1.39 


Similarly 


t  19  jflj.  ;gn£-T*iB  bnB  41. 1.  A  brtfi  >1.1. A  anoir^Bupe  intf 

?.£  9:f  n  )  YE.l.A  nci:ctBJjp9  ectt 


■3 1  o 


A.  12 


Pxz  =  p^a 2ra  Cos  P  {0t_0r*  ”  iPY0i2ra  Sin  P  ( 0^  +  0^ } 

-  pc2y(y-1)  Cos  Q  {ft,+fi  }  +  1pc2y(y-D  Sin  Q  {ft, -ft  } 

T/  ±*  0  jl 

A. 1.40 


A. 1.4.  Matrix  formulation  and  transfer  functions „ 

Equation  A. 1.35*  A.  1.36,  A. 1.39  and  A.1.40  can  be 
written  at  the  in  interface. 


a 


m 


Cos  (e,  +6  }  + 

2  m  t  r 

*  mm 


ia2 

— -  sin  pm  {et  -e_  } 

c2  m  tm  rm 


Y  r  Cos  Q  -nr  )  +  lYmrg  Sin  {J2t  -«r  } 

m  mm  m  mm 


A.  1. 41 


d)  ia  2  a  2 

JH  =  r  sin  P  {0.  +  0  }  -  —  r  Cos  Pm  {0.  -0^  } 

o  a  m  t  r  2  ot  m  t_  r_ 
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Equations  A.1.41,  A. 1.42,  A. 2. 43,  A. 1.44  can  be  expressed 
in  matrix  form: 
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where  the  elements  of  the  matrix  Dm  are 
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(d  )  =  ip  c2y  (y  -1)  Sin  Q 

m  4  3  '  m  'm  m 

(d  )  =  -  p  c2y  (y  -1)  Cos  Q  A. 1.46 

m44  m  m  m  m 

Equation  A. 1.45  can  be  written  for  reflected  and  incident 
waves  which  is  incident  from  the  bottom  of  the  ntJ:1  layer 
(Figure  A . 1 . 5) : 


Figure  A. 1.15.  Direction  of  axis  and  numbering  of  layers. 
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At  the  top  of  the  layer,  dm  =  0  and  equation  A. 1.47  will 
be 
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It  can  be  written  as 
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From  Equations  A. 1.48  and  A. 1.49 
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Matrix  a  =  D  E_1  can  be  calculated  from  equations 
m  m  m  M 


A. 1.48  and  A. 1.51: 
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( am } 2  3  =  i(pmc2rI(ra  Sln  Pm  +  rV  Sln  V 

m  m 
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Repeating  equation  A. 1.52  for  each  interface, 
displacement  at  the  surface  can  be  found  in  terms  of  dis 
placement  at  the  bottom  of  the  layers : 
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From  equation  A. 1.50  and  letting  P  =  P  =  o  at  the 
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Substituting 
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we  can  write 
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In  the  case  of  P  waves  of  unit  dilatational 

amplitude  in  the  layer,  setting  0.  =1  and 

1m 

ft.  =  o  in  equations  A. 1.57,  and  solving: 
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D  =  (J  -J  )(J  -J  )  -  (J  -J  )(J  -J  ) 
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Prom  equation  A. 2. 48,  at  the  top  of  the  layer 
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Taking  the  ratios  of  component  in  equations 
A.1.60,  A.1.61,  A. 1.63  and  A. 1.64,  we  get  Horizontal  and 
Vertical  transfer  functions  at  a  given  frequency: 
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Equations  A. 1.65  and  A. 1.66  are  the  complex  functions  of 
frequency . 

To  normalize  surface  amplitudes  to  unit  total 
amplitude  in  the  incident  wave,  equation  A. 1.65  should  be 
multiplied  by 
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and  equation  A. 1.66  should  be  multiplied  by 
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Using  equations  A.l.63  and  A. 1.64,  transfer  functions 
are  found: 


TV  (w)  =  2(J  -  J  )/y  rft  D 

s  12  22  m 
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TW  (to)  =  2(J  -  J  ) / Y  D 

s  21  1 1  m 
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For  normalization  equation  A. 1.73  should  be  multiplied  by 
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and  equation  A. 1.74  should  be  multiplied  by 
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APPENDIX  2 

LISTING  OF  COMPUTER  PROGRAMS 


This  appendix  contains  the  Fortran  IV  programs 
developed  using  the  theory  in  the  text. 


FORTRAN  IV 

THE  BUTTERWORTH  LOW-PASS 


FILTER  SUBROUTINES 


n  n  n  n  n  n  .  n  n  n  nn 


A.  27 


C 


SUPROIIT  inf  filco 


this  subroutines  calculates  butterwor 


h  low  pass  ultcd  rOEEFiriFMTS 


COMMON /F  IL  /F  (  4  )  /  I  v  T  /  A  T/C  UT  /-°rn 
RF A L  0Oee(4} 

COMDlrX  B  (  a.  )  ,m  (  4  ) 

FRFQ.  T  S  OUT  QFe  F  p  e  q i  |  f  N  r  y 


SR  FORMAT  n  H]  ,  SX  »??HrORBEr  T^o  njj  OF  frfq.,c£#^) 

00  format  (  1  Mi  ,  1  7HCJT  OEF  fPeqijfn^Y  *fs  .  7/ ) 

1  On  format ( F6 . b ) 

101  format ( fb.b } 

1  OR  format  (  RX  » 1  quo  I  G  I  T  I  Z  T  mo  INTFPVAL^b.r/) 

1  0  b  FORMAT ( 1 H 1 »  b  S  w  s - o  l  a  N  f  p  0  L  F  5  FOR  UNIT  LOW  PASS  F  I  LTFR // ) 

1  04  FORMAT  (  1  X  >  FQ  •  A  *  eo  ,  /t  ) 

10S  format  (  1  Mi  ,  i  7HS-DLAMF  ROLFS  FOR  »  FS  •  “ » 1  P.HHZ  LOW  PASS  fjlTFR//) 
10s  FORMAT (  1  HI  , 1 SHCOrFF .  OF  DQLYN.//J 
107  FORMAT  (  1  X  *R  (  IE  ,F-[7.S  )  ) 

1  0  F  FORMAT  (  1H1  »bphf  ILTFP  OOFFFIC  I  F  N  T  S  eqr  LOW  PASS  F°FO  .  ,Fft  .  "*  /  /  ) 

RF A O ( S  » 1 00 1  FReq 
WR  ITF(  S  » 1  on  c R r0 
R F A D ( S  » 1 0 1  )  AT 
T  =  1  .  /  a  T 

W  R  I  T  F  (  G  »  1  0  R  •  T 

F(i  )  =  (_O.QRBO.B.bq->71 
0(R)  =  (_0.Op'5O%_n#7P77) 

S(?)r(-n.Bfi?7)A,n«)1Q) 

S  (  4  )  =  (-0  .  B.QR7  ,  _o  .  o->BO  ) 

WR  I  T  E  (  s  ,  1  0  B  ) 

WR  ITF ( 6  *  10A )  (  B  (  I  )  ,  r  =  1  ,  L  ) 

W  = FR  FQ*  a;  .  R  0  R  •»  QBE 


FRFO.  C ORR  FG T I ON 
X=W*T/R . 

A= { R • 0/T ) *T AN ( X  ) 

V/R  I  T  F  (  6  *  S  ?  )  A 
C 

00  1  1=1*4 

M ( I) =5 ( I ) *A 

1  CONTINUE 

WRITFI6 *  iOS  )  FRFO 
WR ITF (S  » 10A )  ( m (  I  )  , T  =  1  »  4 ) 
C 

1  =  1 

?0  COFF ( T ) =- ( M (  I  )+M(  1+1  )  ) 
COFF ( 1  +  1  )  =M(  I  ) (  I+1J 
1=1  +  2 


1 


A.  28 


IF  (  I  .FO.  F  )  GO  TO  TO 
GO  TO  ?0 
°0  WRITF(6*106) 

WRITF(6*107)  (  I  » G0C  r (  I  )  *  1  =  1  »4) 

C  TPANSFOPO'  IN  TO  Z  -  o  L  A  N  F 
C  ! 

<  =  1 

lio  PFOr?*(;OFF(K+1  ) -P . / ( T  **?') 

AL.P=(4./(T**?)  )  -  (  ?*COFF  (  <  )  )  /T  +  OOFP  {  <  +  1  ) 
rA<=(4./(T**?)  l  +  l^rOFFfk'l  )  /T  +  COFF  (  K  +  l  ) 
FfV'  )=RFD/GA< 

F(F  +  1  ) =  A  L  P / G  A  K 
^  =  vr  +  7 

TF  ( K . FQ . f  )  GO  TO  10O 
GO  TO  110 

ICQ  W  t?  i  t  F  (  A  *  1  f!  p. )  F  R  F  0 

W R ITf(6»107)(I»F(I),t=1»4) 

RETURN 

END 


SURROUTINF  FILTFR(MG) 
rOM^ON/HOLD/D  (TOGO)  ,^d  (  i  toq)  /(tjl/c  (4) 

D T  MFMS I OM  C(1000) 

C  n  )  =PM  ) 

c  (  •>  )  =p  (  ?  )  +D  (  1  )  -F  (  1  )  *r  (  1  ) 

SP ( 1 ) =r ( 1  ) 

Sp ( 7 ) -r  ( ? ) +r ( 1 ) -f ( r ) *s° ( 1 5 

no  in  1  =  f »ng 

r(T)=P(n+P(i-i)-PM)*r(!-!)-F(?)*r(i-?) 
FP (  T  ) =C (  I  >  +  C (  1-1  ) -F  f  ?  )  *FR (  T-l  ) -F ( 4 ) (  I 
1G  COMTINUF 
RFTlJRN 
FNO 


FORTRAN  IV 


THE  BUTTERWORTH  HIGH-PASS 


FILTER  SUBROUTINES 


ro  n  n  n  n  n  n  o  nn  r>.  n  n 


A. 31 


C 


SI.JP ROUT  INF  FILCO 


THIS  SUBROUT  INF  CALrULMcS  RUTTFRMORTh  HIGH  PA  F  S  ^ILTF0  COEFFICIENTS 

rO^ON/F  IL/F(4  )  /  INT/ AT/rijT/cppQ 

REAL  C. O f  p  (  4  ) 

COMPLEX  S(4)*M(4) 


FREO.  IF  T-HF  CUT  OFF  FRFO. 

SR  FORMAT ( 1H1 ,5X * ??HrOPRFCTFn  CUT  0r  CRF0.,F4.^) 

QQ  FORMAT  (  1  Hi  ,  I  7HCUT  OFF  FR  FQUFNCY  *  ’  /  ) 

ion  FORMAT (F 5. 3) 

101  FORMAT <F 8.3) 

10?  FORMAT  (.?  X  ,  1  OHDI  GT  T  I  7  T  NO  I  NT-pRVAL  »  F  p  ,  ^  /  5 

I0*a  FORMAT  (  1  HI  ♦  T8HS-PLANF  PCLCS  FOR  UNIT  LOU  dASf  FILtfR//) 
104FORMAT(1X»FQ.4,Fp.a) 

l  Os  FORMAT  (  1  H]  »  1  7HS-PL.ANF  POLES  FOR  *  F6 . 0  » 1  Q  HHZ  .  H  I  GH  PASS  ft(_TfR//) 

106  FORMAT ( 1 H 1  * l 6HCQPPF .  OF  PQLYN.//) 

1 07  FORMAT (  1  X  •  7 (  T  5  ,  Fl 7  .  S  )  ) 

108  FORMAT ( 1H1  OPHFILTFR  COEFFICIENTS  cOR  HIGH  PASS  FREQ .  ,  F6 . ^ / / ) 

RF A 0(5*100)  FRFO 

W R  I  T F ( 6  »  ° Q )  FRFO 


p  F  A  o  (  s  ,  ]  o  1  )  AT 
T= 1 .n/AT 
WRITE (6, 10?)  T 

S(1  )  =  (-Fi.o?^Q«R.fq?7) 

S  (  ?)  =  (-0.Q?'3O,_n.FQ77) 
F(7)-(-0.FF97, 0,0^0) 
S(4)  =  (-  r'.'58?7,-O.Q?70) 

MR  T  TF I  6  vl 0  7 ) 

WR  ITF ( 6  *  104 )  ( s (  T  )  , I  =  1  *  4 ) 
w  -  F  R  FO  >-  6  •  ?  8  n  1  s  a  7 

FRFO.  CORRFCTION 

X  =  W  *  T  /  ?  . 

A  = (?.0/T )*Tam ( X ) 

WRITE (6 *5?)  A 

DO  1  1=1*4 

MU  )  =  A / S  (  I  ) 

1  CONTINUE 

WRITE (6 *105)  FRFO 

WR I TF ( 6  *  1 04  )  (M (  I)  * !  =  1  *  A  ) 

1  =  1 

?  0  C  O  F  F  (  T  )  =.  -  (  M  (  I  )  +  M  (  t  + 1  )  ) 
COFFC  1  +  1  ) =  M (  I  )*M(  T+l  5 
I  =  I+? 


n  n 
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IF  (  I  .EQ.  *>  )  GO  TO  T  0 

GO  jo  ?0 

U'P  TTF{  ,  1  oo  ) 

V/P  I  T  v  (  ft  ,  i  o 7  )  (  I  »  COcF  (  I  )  ,1  =  1  ,  4  ) 

C 

TRANSFORM  TN  TO  Z -PLANE 
<  =  1 

110  DRP= (  (  T  *  *  ?  ) / ? ) *COF  F ( K  + 1  )-?. 

ALP  =  1  .-(  T/?  )  *rOrF  (<)-•-((  T**?  )  /  A  )  *rOcF.(  <  +  i  ) 

r  ap=  ( T/o  )  *coef  (:<)  +  (  (  t  *  *  ?  )  /  u  )  *rocF  (  k  +  i  )  +  i 

F ( K ) =PFP /CAK 
F  (  <  +  *1  )  =  A  L  D  /  C  A  < 
k'rrk'  +  p 

IF ( K • FQ • S  )  GO  TO  10O 
GO  TO  11 0 

10P  WRTTF(6,10P)  F  R  cO 

WR  ITF (6  » 107  }  (  I  *F ( T  )  » T  =1  *  A ) 

pp  turn 
END 


\ 


. 
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SURROUTINE  F  1 1. T~R  (  NG  ) 

COWION/HOL  n/°(  1000  )  »Gn(  1000) /FIL/r(^) 

D  I  MFNS  TON'  c  (  1  0  00  ) 

c ( n  =r ( i ) 

C ( ? ) =p ( ? )-?*p ( i ) -f ( i ) *r ( i ) 

sr  (  n  =c  n  ) 

SP(?)=r(?)-?*r(i ) -F { i ) *  s  n ( i > 

DO  10  1=0,  MG 

r  (  i )  =p  (  i  )-?*p  (  i  -  i )  4-p  (  i  —  o  )  —  f  ( i  )  (  i-i  )  -f  (  o  )  *r  {  r-o) 

sn( r ) =r ( i ) -?*r ( i-i  ^  +r  ( i  -  ?  )  -f  {  ^  )  *  sp  ( i-i ) -r ( a ) ( 1-7 ) 
10  rOMUN'UF 
RETURN 
FNO 


\ 


I 
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C 

c 

c 


c 

c 

r~ 

V. 

c 

c 

c 

c 

c 


r 


SUPPOUTTNF  FILCH 


this  subroutine  calculates  eu 


T  ER.WOR  TH  BAND  DASS  FILTER 


9 


r  o WON /  F  I  L /  F  I  (  8  )  /  I  N T  /  A  T  /  CUT  /  A  1  »  A  ? 

complex  s  (  a  )  »m  (  a  ) 

real  EOEFF(q) 

i'  1  ,  A  ?  APE  LOW  A  N  O  HIGH  rUT-Q"F  FREQUENCIES 
1/AT=T  DIGITIZING  INTERVAL 

S..  IS  S-PLANe  POLrS  OE  LO’,;  PASS  R'  JT  T r  R  WORTH  FILTER 
FT  (  P  )  .  .  are  EJLTER  rncFHr  I -NTS 

COFEF.  OF  POLYN.  . . •  ARP  C  0  c  E  p ! G I F  N  T  OE  POLYNOMIAL  OF  TRANSFER 
F 1 1 N  r  T ION  IN  e-plane 

ICO  FORMAT  {  ?  e  q  #  -3  ) 

TOl  format  (  I  H I  »  t  o  h  Z  -  D  L  A  N  e  dri_fS  cor  LOW  PASS  FJLTfp//) 

1  Or  format  (  IX »F7.4 Jf7, A ) 

10-3  FORMAT  (  1  HI  »SX  >  RTHCORREg  T  eo  CUT  OFP  FREQ.//) 

104  FORMAT { 1  OX i ?F 1 n . ? ) 

1  0  E  F  o  R  M  A  T  (  1  H I  » 1  X  ,  1  s  H  S  p  L  A  A.;  F  POLES//) 

j  0  6  FORMAT ( I  X  » 1  (  I S  , F 1 R . 4  -  p 1 ? . 4 )  5 

107  FORMAT { 1H] *SX » 1 6uCOFFF.  oe  polvm.// ) 

TOP.  EORMAT  (  1UI  >  EX  ,  ?7Mf  iltpp  GO-EP.  ear,  .  .  .pp  .9  ,  -|  X  ,ES  .?// )  ’ 

1  0  0  E  o  p  a  T  (  1  X  ,  ?  (  IS»1X*Ei?.a)  ) 

S  (  1  )  =  {-  0,op?o,o,^.q9-7) 

S(?)=(-0.p?^p»-0.7q?7) 

S(4)r(-nOQ^,^-tOMO 


c 

R  p  A  D ( e j 1  00)  A  1  ,  A  ? 

UP T  T  F ( A  »  1  0  1  ) 

VR I T  E ( A  *  1 o  ?  )  ( s (  I  )  ,  T  =  ^  '  4  ) 

c 

T= 1 .a/at 

V,  1  =  A  1  *  6  •  B  p  E  i  O  5  o 

W  ?  =  A  ?  *  C  .  7  P  7  1  ft  S  7 

c 

X  =  VM  *T/?  . 

Y=W?*T/?. 

C 

A  =  (  ?  .  /  T  )  *  T  A  N  (  X  ) 

R= { ? . /T  } *T AN ( Y  ) 

UR  I  T  F ( S  » 1 0  e  ) 

V/R  ITF(  6  »  ln4  )  A^e 
A  A  =  B  -  A 
P  R  =  A  *  R. 

C 

on  7  0  I  - 1  *4 

v  (  j  )  =  (  s  (  I  )  *  A  A  )  /  ?  .  -  G  FOR  T  (  (  (  4  A  *  c  (  I  )  /  7  .  )  *  *  7  )  -q  r  ) 
v (  j  +4 )  =  ( S  (  I  ) * A  A ) / 7 . +G SOR  T (  (  ( A  A  *  S (  I  ) / ? .  ) ** 7 ) -p r ) 
70  CONTINUE 


' 


c 


'<'R  !  T  F ( F  » 1 05 ) 

WR  I  T  F  (  F ,  1 06  }  (  (  I  ,  m  (  I  )  )  ,  T  =  1  ,  q  ) 


I  =  ! 

?n  rOF^F  (  n  =-(M(  I  )+\H  ui  )  ) 

COFFF  (  T  +  1  )  =  (  I  )  (  I  +  1  ) 

I  =  I  +  ? 

IF(-I.FQ.o)  r,C  TO  ^>0 
GO  TO  ?0 
10  V/R  ITF-(F>  1  07) 

WR  I  T  F  (  A  ,  ]  0  Q  )  (  (  I  ,  rri  r  c  c  (  j  )  )  .  t  -  i  ,  p  ) 

Y-  1 

fo  iL°=  (  GOf^f  ( '<+1  )  *T  )  /?  .  -rocrFc  {  ^  )  +>  .  /r 
RFD  =  rOP'c  F  (  K+1  )*T-4./T 

C*A,<=(rOfr“F(l/  +  1  )*T)  /  ?.+rOrrFtr(K)+?./T 

r  T  <  <}  =  RFP/C AK 
F I  ( <  + 1 )  =  a  i.  p  /  r  a  f 

<  =  F  +  ? 

I  F  (  F  •  GO  •  °  )  r-0  TO  A  O 

GO  TO  FO 

A  0  WR I T  f ( F  » I  0 3 )  A  1  )  M 

V/R  I  T  ^  (  F  *  1  °  O  )  (  K  *  r  I  (  <  )  ,  <  =  1  •  °  ) 

RF  TURN 
r  N  0 


* 


t 


SUBROUTINE  PEVERS(N) 
COMMON /TUM /SB ( 6S°n  ) 
NN=N/? 

DO  TO  I=1,NN 
J  =  N-I 

TEMP=SR ( T ) 

SR ( I )=SR ( J+l ) 

SR ( J  +  l ) =  T  F  M  P 
in  CONTINUE 
RETURN 
END 


SUBROUTINE  SKIP(N) 
DO  S  K-\  »N 
S  RFAO(A) 

RETURN 

END 


SUPPOUT  I  ME  PEf-'OVF  (  MM  » NG*MNV  »NZ  ) 
COMMON/HOLD/P  (  6  ) 

DO  1  L= 1  »NN 
GtlM  =  o#n 

DO  ?  1=1 » MG 
?  5UV=SUM+P(L» I ) 

A V  =  SUM /FLOAT ( MNN ) 

WRITE (6  ) 

V/R  T  T  E  (  6  »  A  )  A  V 
DO  E  K  =  1  ,MMM 
P(L»K)=D (L»K)-AV 
E  C  O  M  T  [  M  U  r 
1  CONTINUE 

^  FORM  AT  (  I  OX  ,  T^HAVEP  AGF  VALIJF/) 

A  FORMAT ( 10X  *  F20 . ? ) 

RETURN 

END 


' 


SUP  ROUT  IMF  FILTF9(L*NC?) 

HTMcmsTON  r  (**  )  ,D(  *>  )  *P  (“*  ) 

COMMON /HOLD  /P  (  6  )  /tijm/Fd  (  )  /p  JL/F  I  (  R  ) 

MOn^  (  INOFX  )  =  I  MOF'X-  (  (  I  N  R  c  X  - 1  )  /^  )  #7 
rm=°(Lti] 

c ( ? ) =D ( l  *? ) -ft ( i  )*r  (  n 

C  (  ? ) =p ( L  ) -P ( L  ♦  1  ) -F  I  m  *c ( 2 ) -F I ( ? ) *C (  1  ) 

D(1  )=C(  1  ) 

D(  ? ) =c ( ? )-F I ( p )*n( 1  ) 

D('a)=C(^)-C(l)-pI(^)*n(?)-crT  (41*0(1  5 

f  (  n  =  o  ( i ) 

F(?)=D(?)-FT(F)*cr{1  ) 

p(7)=n(p)-r)(i  )-Fi(s)-^p(9)-Fi(6)*crii  > 

SR ( 1  )  = F (  1  ) 

SO ( ? ) =C ( 2 ) -F I ( 7 ) *sp ( 1 ) 

SR ( ^ ) =F ( R ) -F (  1  )  -F  T  { 7  > *SR ( ? ) -F I { R ) * SR ( 1  ) 
no  10  I =4* MG 
I  l  =  W0Dq (  I  ) 

imil=mods  ( i-n 

IM?L=wODP ( T-? ) 

C (  IL ) =p ( L  *  I) -p ( L  * T-? I -FI  U  ) *C ( I  MIL) -FI ( ? ) *r ( TM?L ) 
0  (  IL  )  =C  (  IL  )-C  (  I M ^ L  )-c  T  (  3  )  *0  (  TM1L  )~c  I  (4  )  *r>(  IM?L  ) 

F(  IL  I  =n  {  IL)  -0(  IM?L  )  -F  I  (  5  )*F  (  JN’lL  )-f  I  (  S  )  *  F  (  I  w?L) 

SP  (  I  )  =F  (  IL  )-F  (  I  VI ?L  )  --  I  (7  )  *SP  I  I-l  )-=•  I  (P  )  *SR  (  I-?  ) 

10  CONTINUF 
RF  TURN 
END 


. 


S M D P 0 U T  IMF  PFAn(lPFC)' 
COMMON /HOLD /D (6*6^00) 
COMMON /  T>  JM/  I  A  (  1200  ) 

I  N  D  F  X  =  0 

DO  10  J=1 »LPFC 
R  F  A  D ( 4 )  IA 
00  2 n  Ln  =  l *6 

I N D f X  =  ? 0 0 * ( J-l  ) 

DO  20  lRP=LP»1200»6 
I  N DF X  =  T  NDFX  + 1 
?o  P  ( t  f  ,  IMDFX  )  =- (  I  A  ( l.RO  )  ) 
10  CONTINUF 
RF  TURN 
FND 


A.^2- 


SUP-ROUT  I  ME  DDl_OT  (  MM  ,wnm  » 3  U  c  » T  p  I  L  *  T  R  p  p  ) 

COMMON /.HOLD /P  (  6  » )  /TUM/SP  (  5,^00  )  /  p  IL/p  [  (  Q-  ) 

COMMON  /  CONST  /T  IN, UN  I  T  *  AMP  ,  AMPMAX  ,CSFP  *  SHOT  ,  TCOR  *  p,  /  INT/T/CUT/A1  *A  2 
RF C  =  FLOAT (  I  RFC ) 

XSCALF=AMP/ AMPMAX 
Y S C  A L F=  T  I N /UN  T  T  i 

CALL  PLOT(n.o,n.n,'5) 

CALL  SVMPQL  (noE0»-o.^0,P.?0«i,7HR^COR0  A  ,0.0,]?) 

CALL  NUWPFP  (  i  .  7  ,-o  .  ^  ,o.?o  ,RFC»n.  "  ,  -1  ) 

CALL  MUMPER  (?.5*-0.3O,0.;>0*rW0.0',-!  ) 

CALL  5  Y  M  P  0  L (D. 5 0,-1 .  oo  ,  "  .  2  o  ,  ahSHOT., O.o,4) 

CALL  NUYRER ( 1 .5,-1. 00, O.?o, SHOT, 0.0,1) 

CALL  SY^POL ( 0. 5  0 ,~1  .6° >  ^ . 15 , 10HT  p'P  CORR . ,0 .0 , 10 ) 

CALL  N'JMPCR  (  1  .o  ,-l  .fto  ,o.  1  0  ,  TfOR  ,  °  •  p  ,  ?  5 

CALL  SYMPQL ( 0.50  1 0,0.1s , 10HDIGIT . INT. ,0.n , io ) 

CALL  MUMPER ( 2 .0 •  -  ? . i o  ,0.10  ,  T  ,0.0,4) 

CALL  t'YMDOL  (0.53,-2.60,0.1  3  ,  PH  ARC  LTNP,n.n,a) 

CALL  SYMROL  (0.50, .  i  0,0.1c  ,omPpojpcT  °,o.r',o) 

CALL  SYMBOL  (0.50,-**  .50.0.15 » 1 PHr  ARL  Y  R  I  SF  ,0  .r> ,  io  ) 

IF ( IF IL.FO. 1 )  GO  TO  ] 0 

CALL  SY^POL  (  0  .  =10  ,-4 . 1  0  «  n  .  1  5  •  10HR  TLTFR  RV.,P.on,in) 

CALL  SYMBOL  (0.50  , -4 .60 ,0.  i  R  ,‘15H  HZ, ,0.00, IB) 

CALL  NUMRFR (0.50,-4.60,0.1 b ,ai ,o.p ,2 ) 

CALL  NUvrtp ( 1 . 70 ,-4 .6° , n . 1 p • A? • 0 .o ,7 } 

10  CONTINUE 

TUNIT=IFIX  (LIMIT  ) 

DO  1  J=1 *MN 

DO  ?  1  =  1 ,NNN* I  UN  IT 

CALL  SYMROL ( ( { FLOAT ( I -1 ) *YSC  ALE  >  +  3 . 0  ) ,0.0,0.12,3,0.0,-!} 

CALL  PLOT  ({(  FLOAT  (  I -1  )  *YSCALp  )  +  ? .  0  }  ,0  {  j  ,  j  )  *x SCALP  ,  0  ) 

JSTOP=I+IUNIT 

IF  (  J  S  T  0  D  •  G  T  .  N  M  M  )  C-o  TO  4 

DO  3  I  I = I  ,  JRTOP 

CALL  PLOT  (  (  (  FLOAT  {  I  I  -  1  )  *YSP  A  Lp  )  +  '>  •  0  )  ,P  (  J,  I  T  )  *XSC  ALC  *?  5 

3  CONTINUE 
?  CONTINUE 

4  KKsJSTOP— IUNIT 

CALI  PLOT  (  (  (PLOAT  (’<’<-1  )*YSCALc)+3.0)  ,D  (  J  ,<K  )  ^v^^alp  ^  I 
DO  5  <  =  «,N\'N 

CALL  PLOT  (  (.(FLOAT  (  K  -  1  )  #YSr  A  LF  )  +3 , 0  }  ,PfJ,K)*X5rALF,?) 

6  CONTI. NUP 

CALL  Pt.OT(0.o,CSPO,-3) 

1  CONTI NUP 

CALL  PLOT (0.0 ,-0.3  , -3  ) 

DO  6  1  =  1  , NNN , I  UN  I T 

CALL  SYMPOL (((FLOAT(l-l)*YSCALF)+3.o),0.n,o.^o,iB,o.o,-l) 

A=FLOAT ( 1-1 ) /117. 

CALL  NU^rfR  (  (  (  FLOAT  (  I  -1  )  *vsfALF  )  .0  )  ,~0.3o  ,0 .00  ,,a  ,0,0 ,0  ) 

6  CONTINUE 

CALL  PLOT(0.o,4*CSFP,-3) 

RETURN 

END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


BUTTERWORTH  BAND  PAS?  FI|_TEp . . . 

THIS  PROGRAM  CALCULATES  FILTER  COEFFICIENTS  AND  F I L  T  F  P  ANY 
DIGITIZED  IN  PUT  .  CONSTANTS  TO  Rc  rear  ARF  ON  COMMENT  CARD* 
IN  EACH  SUBROUTINE. 


MD.  .  .M!J'-’«ER  DF  BLOCKS  TO  RE  SKIPdpq  ON  TAPE 
NF  •  •  .NLJvaER  R  L  0  C  F  e  TO  Rear 

NF,..mUM°er  oe  q  rrop^.S  j  r,  qp  PROCFERpn 

nG.*.NUmqer  of  w o  r  p  s  per  oh am vfl 

max?  MUV  a  l  L.  O  W  a  R  L  c  MG  IE  ftqpq 

MN.  .  .NU'-nFP  OE  CHANNELS  TO  »F  FILTERED 

NNN.  .NIJ^RFR  of  WORDS »  exCLDS  IVE  of  LEADING  AND  TERMINATING  zeros 
HZ...NUHEFD  qp  Z  e  p  O  S  after  BLOCK  oe  WORDS 

!OL0T  =  l  dLOT  THE  DATA  RpFORF  pILTeo  -  . 

=  0  DO  NOT  PLOT  THE  rata  q  E  P  0  R  c  FILTER 
IW  =  1  V! D  I  T  E  FILTERS  OUT  PUT 

=  0  DO  NOT  ’'.'RITE  EIlTeRpo  pm  it  OUT 
<L  I  ST  =  0  DO  NOT  write  TV  -5 1 1  r  data 
=  1  WRITE  THE  IN  OiJT  DATA 


COMMON /HOLD/0  (6  »  f  ^  0  °  )  /  T,|VVSn  (A^^o)  /  F  I  L  /  e  T  (°  5 

COWON /CONST  /T  I  N  « ■  IN  T  t  ,  a  mo  ,  \  vdvax  ,  r  -po  ,  SuOT  JSOR  ♦  °  /  I  N T /  T 

REAL* A  FUF(SOO) 

KK=A*500 

SFT  THE  WORK  AREA.  FOP  PLOTTING 
CALL  PLOTS  (  a U E  (  1  )  *  V’K) 

CALL  PL0T(S^en,?7.^-^) 

READ  CONSTANTS 


RPAn(S,100)  MO  »NF  ,NF  *  mn  »  I  PLOT  •  I  W*KL  T  S  t  ,  NNN  «N7 
RFAD(S*10?)  Gain 
WRITE (  Fm  19° ) 

WR  ITF(6»l°f>)  GAIN 

Read  (  S  *  1  05  j  T  I  N  »UN  I  T  »  AMP  ,  AM°maX 
WRITE (s , 104 ) 

V/P  ITF  (6  »  1°5  )  TIN  *UN  I  T  ,  AMP  »  A  V o *.*  A  X  »  C S E P 
N  G  =  M  N  N  +  m  Z 
WR ITE (*  ,  9 Op  ) 

WR  ITF ( 6  ,  ROo  )  ( NO  ,  me  ,NE  ,NG  »NN « I  PLOT  ? IW,KL T ST , NNN ,N7 ) 

CALCULATF  FILTER  COpepigIpNTS 

CALL  F I LCO 
C 


- 


■*- 


' 


n  pi  non  n  n  n  no-  n  non  on 


A.  <45 


I F ( NO • rQ • o )  GO  TO  l 5 
CALL  SK I ° ( NO ) 

IS  CONTINUE 
C 

DO  ?  50  IL  =  1*NF 

C  PF.AO  RFCOP  0  PARfi^FT^P 

C 

C  SHOT  TF  FH^T  NU^pfR 
C  TCOR  IF  T  t  k<  c  C0PR  FCT  ION 
C  R  IS  A p C  LIN-  NU-VRFR 

C 

REAP)  (  5  »  1  OF  1  SH0T»TC0R*P 


IRFO=  (  ND  +  NJE+  (  I  L  -  1  )  *  0  P  )  /  F  p 
I F (  I ° r c . c Q . ° )  lDrr  =  i 
TF(MF.LT.PO)  IR-^=(\!0/?P)+T 

PUT  ZEROS 

DO  1  L  =  1  » N N 
DO  ?  J= 1 » NG 
P  (  L  *  J )  =  0 . 0 
p  COMT I MUF 
1  CONTINUE 

CALL  R F  A  0 ( M F 1 

WR  I  Tp ( 6 > 004 ) 

W R  I  TP ( 5  * ?CA ) 

CALL  RPHOVF ( N N »NG  »NNN »NZ ) 

I  F  ( V' L  I  S T  .  F 0 . 0  )  Gn  TO  ?? 

00  40  LF=] *nn 
WR I TE ( 6  *  1 007 )  LF 

40  WR I TF ( 6  ♦  I  0^5 )  ( D ( Lc  «  Lr  )  *Lr=l  »  MG ) 
22  IF (  I  PLOT  .  FO • r )  GO  TO  * 


I F  I  L  =  1 

CALL  OPLOT  (  MM  *  MMM  »D,fF*IFTL»I  7rr  ) 

VP  I  TP (6 . 1 ) 

VR I T  F ( 6 ♦ 1 0  o  ^ ) 

F  CONTINUE 
DO  ?5  L  =  1 » MN 

FILTFR (((((((((  (  ( ( (  (  (  (  (  (  (  (  (  (  (  (  U((  f  f 
CALL  EILTFR(L*NG> 

CALL  RPVERS(MG) 

C 


A  .46 


DO  6  J=1 *NG 
P(L*J)=SR( J) 

6  CO^TINUF 

CALL  FTLTFP(LvNG) 

CALL  RFVFRS(NG) 

00  1?  J=!*NG 
P  (  L  *  J )  =  S  n  (  J  ) 

P(L»J)=D(L»J) /GA  I  V 
13  CONTINUE 

IF(IV#cq.o)  GO  TO  33 

c 

UR riF(3 , 1007)  L 
VPITF(6*201) 

WR I T  F ( 6  «  200  )  ( R ( L  *  J  )  *  J  =  ! »  MMN ) 

23  CONTI  NU  F 

C  (((((((((((  M  (((((((((  f  (((((((((((  (  ( 

C 

IF  Il  =  o 

CALL  DP  LOT  (  NM  ,  NNN  ,  Rl  IF  ,  I  F  T  L  •  T  RCC  ) 

WRITF(6*10^Oj 
WR ITF ( 6  » 10  00) 

C 

230  CONTINUE 
C 

100  FORMAT (RI6) 

10?  FORMAT {^33.1  ) 

103  format  (  TF.g  .  3  ) 

‘104  FORMAT  (  1  HI  ,  1  X  »  4  1  H - TIM - "NIT - A  mp —  A'-^Miy 

103  FORMAT ( f  f  8 . 1  ) 

1  q6  format ( IX  ,fi B  .  1  ) 

IRQ  FORMAT ( 1 H 1 *4HGA  IN// > 

2  00  FORMAT (  1 0( IX »F1 3.3  )  ) 

201  FORMAT ( 30X O ,MF I NAL  OUT  RUT//) 

2np  fORMat  M  Ml  ,iy,7?H - NQ - NF - NF - NG - NN- 

. — KL  I  r  r - NNN - N? - /) 

2  03  FORMAT  UX  *1  3  n  X  «  !6  )  ) 

204  FORMAT  (  1  X  OH  /  /  ) 

1000  FORMAT  (  1  X  ,?7HPL0TT  TNG  HA  8  RFF.N  COMPLFtfr) 

1 006  FORMAT ( 1 X »: PF7 . A ) 

1007  FORMAT (iohi  CHANNfl  , !4) 

REWIND . 4 

CALL  PLOT(n.n,o,n,ooq) 

3  TOP 
FND 


% 


03 ep/ ) 


-iplot- 


'  ■  . 


FORTRAN  IV 


FILTER  IMPULSE  AND 
AMPLITUDE  RESPONSE  PROGRAM 
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c  This  program  CALCULATFS.  IMPHLSO  AND  a  m  p  l  I  T  t  j  d  F  RFSPONSF  OF 

C'  FILTER 
C 

C  NOPTS . . • M  ?  J  * 4  =  F  R  OF  POINTS  i 

C  NP • . . NUMRF R  OF  POINTS  TO  PLOT 
C 

DIMENSION  A  Y  (  ,fX  (  p^n)  ,  ''"Y  (  PAn)  ,qV  (  p?  )  ,  ov  (  p?)  ,  Rt. 

COMMON /HOLD /P ( 1000)  , cr ( 1 ono )  / T  N T / A T /CON ST /NOPTS 
COMMON /DAT /X (  4  A  0  )  , Y ( 4 « 0 ) 

C 

c 

qo  FOPMAT  (  ipil  p  OX  » ISHlMPlJLSo  R c SPON Sr  / /  ) 

1 Oo  FORMAT ( 0  I A ,co . 1  ) 

10?  FORMAT  (  l  HI  ,  1  9HRFAL  AND  I  ma,G  PARTS//) 

1  0  A  FORM  A  T  (  1  X  *  A  {  fSP.XGM.?)  5 
1O4  FORMAT  (  1H1  ,PHF  I  L.TFP  1//) 

1  OS  FORMAT  (  1  hi  » .QHF  I  LTFP  ?//) 

TOG  FORMAT (  1 X  >4 (  I F , 1 X ♦ 9 . p , I  X  »  FI R . T )  ) 

107  FORMAT  (  1H1  POX,  OHAMPL  I  TURF//  ) 

1  0  o  F  o  P  M  a  T  (  1  X  9  r-  (  I  F  »  1  X  »  F  1  F  .  I  )  ) 

1  0  Q  f 0 R M  a  T  (  1  H 1  ,  F  X  9  4  HG  A  J  N  /  /  ) 

1  1 0  FORMAT ( FX ,F70 . f ) 

C 

RFAD(5»ln0)  NOPTS,  njp,  AT 
C 

K'K=?000*4 

GALL'  p  L  0  T  S  (  D I.  i  F  (  1  )  ,  m  V  ) 

c 

C 

CALL  filch 

CALL  Z f P O ( N O p T  S  *  D ) 

c 

P (  ?  4 1  )  = 1  .0 

GALL.  FTLTFR(NODTS) 

c 

WR I T  F ( A  « 1 ° 4 ) 

MR  I TF ( 6  » 1 0? )  (  (  I «  SR  (  I  )  )  » I =]  • NOPTS ) 

C 

CALL.  RFVERS  (  NOPTS  ,  fr  ) 

CALL  MOV f ( NOP  T  S  »  Sp  «  P  ) 

C 

GALL  F I LTFR ( NOPTS ) 

WR T  T  F ( A  i  1 0  5  ) 

VR I T  F ( A  *  1 0  0 )  (  { T  *  S  R (  I )  ) » I = 1  *  N0D  T S ) 

C 

CALL  RFVFPS ( N0DTS , SR ) 

C 

WRITF(6  *°R) 

V'R  I  TF  (  A  i  1  03  )  (  (  I  »  SR  (  I  )  )  *  T  =  1.  » NOPTS  ) 

C 

C  PLOT  IMPULSE  PFSPONSF 

C 


RUT T FRWORTH 


I F  (  ■?  r\  n  0  ) 


- 


DO  4^0  1=1  » MOPTS 
C  Y  (  I  )  =  S  P  (  I  ) 

rx (  T  )  =  F  LO A  T (  I-]  ) / (  A  T  —  1 ) 

4  F  0  CONTINUE 
C 

CALL  WDLOT  (  CXiCYi  VOD  t  c  ,nijF  ) 

C 

C  FOUR  I FR  TRANSFORM 

C 

CALL  MnVF  (  MODI'; ,  fp  ,X  ) 

CALL  ZFPO ( MOOTS , v ) 

c 

K  =  N!OPTS/?  +  l 
L=NOPTS/2-l 
C 

DO  IF  I  =  1  *L 
X  (  I  +  K  )  =  X  ( K  -  I  ) 

IF  CONTINUF 
C 

CALL  FOUR  ID 
CALL  SORT  ID 
C 

WR  I  T  F  (  f>  »  1  o  ?  ) 

WR  I  TF ( A , 1  OF  )  (  ( JT  »  X ( JT  )  * Y ( JT )  )  «  JT  =  1  ? N0D  T f  ) 

C 

C 

no  *0  J T  =  1  ,  MOPTS 

AY ( JT ) =  S  O  R  T ( X ( JT ) **?  +  Y ( JT )**?) 

3  0  C  0  N  T I N  U  r 

WR I T  F ( 6  » 1 0 7 ) 

WR  T  T  F ( 6  » 1 )  (  (  JT  » A  Y ( J  T  )  )  ,JT=!«N0PTF) 

C 

CALL  max (NODTS* ay ,XmaX  ) 

WR I T  F ( fi  ,  1 r  O ) 

W  P  T  T  F  (  6 , 1 1  0  }  x  m  A  v 

C 

c 

DO  170  1=1,100 
py ( i ) = ay ( I ) 

SFC= AT /FLOAT ( NOPTR ) 

PX  (’n=SFC*FLOAT  (  1-1  ) 

170  CONTINUE 

c 

CALL  AdLAT(rX»pY»pUf»M0) 

c 

CALL  PLOT ( ^.0,0.0 ,ono ) 

C 

STOP 

FND 


SlIP  ROUT  I  NF  ZPRO(LX.X) 
D  I  N*  P  N  S  TOM  X  (  1  0  ) 

I  F  {  L X  )  SO  *SO , 1.0 
10  DO  ?0  I  =  1  ♦  L X 
X(  I  )=0.0 
?0  CONTINUE 
so  RETURN 
END 
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SIJRROIJTINF  MOVF(LX»r  *  0  ) 
0  I MFNS ION  C ( 1 0  )  »  H ( 1 0 ) 

no  io  i  =  i  ,lx 

0 (  I ) =C (  I  ) 

10  CONTINUE 
RETURN 
END 


. 
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SllPPniJTlWF  RcrVrPS  (M  »X  ) 
niVFNSTOM  X(MPPTS1 

f>'M  -  M  /  ? 

DO  ]0  !=1»NM 
J  =  N-I 
TFmd=X ( 1 ) 

X ( I ) =X ( J  +  1 ) 

X ( J+1 )=Tfvp 
10  rON'TINUF 
RETURN 
FMH 


SUBROUTINE  FILTFR(NG) 

D  I  M  F  N  S  I  O  N  cm.  D(8 )  ( b  ) 

-  COMMON / HOL 0 / P (  1  0 n o )  , S R ( 1 0 0 " ) / r I L / F I  ( 8  5 
MOr)MINn*rX)  =  INpFX-((TNOFX-l)/‘a)*‘3 

c ( i ) =p  (  n 

C ( ? ) =P ( 7  )-F  I  ( 1  )  *c ( n 

r  {  b  )  =F>  (  -a  )  -P  (  1  )-FKi)#r(?)-'KN*r(l  ) 

DUKll  ) 

D  (  ?)=C(  ?  ) —  F I  ( 8  I  *  o (  1  ) 

o ( 7 ) =c ( b ) -c ( 1 ) -F I { 8 ) *n ( ? ) -c I ( 4) (  1  ) 

F(i)=on  ) 

F ( ? ) =0 ( 9  )-f T  ( 8  ) *f (  1  ) 

F(*»)=:0M)-nM  )-Fi(U*F(?)-q(A)^(i  ) 

SB (  1  )  =F  (  1  ) 

SB ( 7 ) =F ( 2 ) -F I ( 7 ) *SB ( 1 ) 

SB ( b ) =F ( b ) -F ( l  ) -f I  ( 7 ) *$B ( ? ) -F  r ( 8 ) *SR ( 1  ) 

DO  10  I =4 9  NO 
I l  =V0D8 (  ! ) 

IM1L=V0D8 ( 1-1 ) 

I M  ?  |  =^ODB (  I-? ) 

0  (  IL  )  =p  (  I  )  -p  (  I-?  )-cI  (  1  )  *C  (  I  "1  L  )  -F  T  (  ?  )*C  (  IM?t.  ) 

D  (  !  I  )  -r  (  7  L  1  - r  (  T  M  ?  L  )  -  C  T  (  Q  }  * d  (  I  M  l  L  )  -  r  T  (  4  )  * D  (  T  | } 

f  (  I L  )  =  n  (  T  L  )  -  n  (  IN  n  }  -  F  T  (  B  )  *  F  (  7  V.  1  L  )  _  tr  r  (  f,  )  *  f  /  j  v  o  (..  ) 

SR ( I  ) =F (  IL ) -F  C  I M 9 L ) -c I  (71 *SD ( T-l  )-F  T ( 8 ) *SR (  I-?  ) 

0  CONTI  MU f 
RFTURN 
FND 
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SUPROUT  I  Nc  V1 D L 0 T  (  rx  ,rY  ,  N  *PLF  ) 
pTVFNSIO'm  CX(  1  o 0 n  )  . r y  (  i^o ri) 
rn'<vQM/rUT/ a  1  ,i?/rOMr-  T/NO°TF /  TNT  /  AT 

CALL  DL^T  (n#n,n.n,v, 

CALL  PLOT ( o#n ,6 • ^  *  ? )  i. 

CALL  PLOT( lo.n ,s.S ,9  5 
CALL  PLOT  ( 

CALL  PLOT  (r',0O,°t?) 

CALL  SY'-^OL  (  1 . 0  *  A  .0  >  °  .  1  ?  >4RU  P-'PULSR  RFSPONSF  OF  R!  J  T  TFRVOR 
.ASF  FILT-P  ,^.o  ,ia  ) 

CALL  SY'-*nOL  (  S  .  0  .  *  ,  o  .  t  o  .  i  -*h  T  N  D,iT  d  a  o  A  vct  rp  c  ,  o  .  o  ,  1 7  ) 

CALL  SYv°OL  (  S  .  n  ,/i .  A  ,  o  .  ,ig  v  =  ]  ,  ^  ^  j 
CALL  SVMPOL  (  f>  •  ^  •  a.  .  A  .0,0.0  «pt~|  .  o  .  o  ,  «  ) 

C  A  l  L  NU ^ppp ( 6  .  Q  ,  A  .  6  ♦  o  .  *  q  ,  A  1  •  0 . 0 , 1  ) 

CALL  SYVRQL(F.o,.l.^  o.os  ,RhT=  ,  o  .  o  ,  o  ) 

T  =  1 » /AT 

CALL  NU'^po  (  5  .  ?  ,  l,  .  o  ,  «  .  op  ,  T  ,  "  .  o  ,  a ) 

CALL  SY'^POL  (A  .  F  ,/t  .q  ,o  .  aa  ,  °HHCr=  «.  o  .  o  ,  q  ) 

FALL  M  l !  p  F  P  (6.P»A. 0,0.0  0  ,^,0.0,1  ) 
p  =  F  L  0  A  T  (  M  o  d  T  S  ) 

CALL  SY^nOL(ci.A.i,^,r>,^Q.^''!=  .  ^  .  0  ,  c  ) 

FALL  W  ) va  PR ( f . ? , 4 . o , o  # o q  ,  o , o „ o , i  } 

CALL  PLOT  (  O.S  ,  o  ,  q  ,  p-a  ) 

CALL  SCALF(FX,9.°/'!  »  1  »  1  °.n) 

CALL  SC  a  l  F  t  F  Y  *  5 . f  ,  N  ,  1  ,  1 n . o ) 

I  L  =  N  +  ? 

TY=m+1  • 

FALL  AX  I  S  (  0  .  o  ,  0  .  O  ,  OH  *  VDL  T  T'  »r>F  ,Q,c#p,oo.o,fY(TV’)  *  FY  (  T  L  )  *  1  0  .  " 
FALL  AXIC-  (0#o,rt.o,^MrcrOMO»f>«o,o,/'.o,ry(T<)  ,rv(  nj  .  i  *  ,  o  ) 
CALL  L  T  N  P ( c  X  *  r  Y  *  V , ! *  o  ,  o  ) 

CALL  PLOT { ° . p  ^  lO.o.-o  ) 

\>!R  I  T  f  (  s  ^  ?  ) 

VP  I T  E ( 6  * ?  ) 

p  FORMAT  (  1  X  »  ?7~'DL0TT  INC  HAS  RFCA'  C0VP[.FTF05 
PFTURN 
FA1 0 


BAND 


n  n  n 
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‘'HP  POUT  IMF  *dL0T(PX,d.Y  *ouf  ,.\|P  ) 

c . . . 

c 

DTMFNSTON  RX(  10?  )  *c»Y(  10?  ) 

CO^MOM/rUT//'  T  >  A? 

ORAW  THC  OUT  LINp  Oc  THF  PAT,f 

CALL  PLOT  (  o  .0  ,  0 . 0  « *5  ) 
call  DL^T{0,n,6.^,?l 
CALL  PL°T  ( 1  n,0i6.!i  u) 

CALL  PLOT(lo.n*o.o*?) 

CAL!  PLOT(0. 0*0.0,?) 

CALL  SYMROL  (  1  *0  *fi  .0,0  .  i?  ,a,pha.vpl  TT!  iOF  RFSPQNSF  OF  BUTTERWOR’TH  RANn 
.  PASS  F  1  L  T  P  °  FOP  -  HZ  *0.0,  ft  «)) 

CALL  NU^RFR  *  o  .  o  ,  i  ) 

CALL  M'JVia.trp  (7.B,ft.0*O.1?*A?.Q.''*l  ) 

c 

CALL  dL0T(0.p,o#^,7P) 

C 

CALL  S  C  A  L.  p  (  c  X  »  o  •  ^  »  v  d  •  1  *40.  o  ) 

CALL  S  C  A  L  F  (  p  Y  ,  S  ,  ^  P  »  1  »  4  n  ,  °  ) 

C 

CALL  AX  I  S  (  0 , 0 . 0  •  0  ,  qh a  ,JIPL  I  TUpc  ,  o  .  p  .  p  ,no  .  r  ,  pv  (  i  )  ,  p  v  n  np  )  .  ^  .  r  ) 

TALL  AXIS(/',0,n.^,OHrD-Q'!FNOv.P.P.'''.r'r'.r'.PX(iPl).°Y(inp),60.n) 
oall  L  I  Nc  (  RX  .  °Y  ,A’P  .  1  ,  ^  ,  0  ) 

VP  I  TP ( ft  ♦  ? ) 

WR I T  F ( 6  » ? ) 

?  FORMAT  (  1  X  »  ?7HDL0TT  r  VC-  HA  S  PFFN  CCVPlftpo) 

RF  TURN 
FND 


' 

n  n  r\  n  n  nnnnnnnnn  n  n  n 


SUDROUT  I  ,N!F  FOUR  in 
CNF  D  I  Y  c  r  1 F  1 0 0  A L  roiJP  f  cl? 


A. 56 


C ^VMOM  /m,vc,T  /NOD  T  <; 

COMMON  /  OAT /X  {  4S0  )  <  v  ( ) 

TNTFGFR  J  v<  ,  v  ,vp  ,J!i  J?  >J°  >  JA  ♦  J-  *  JT 
RfAL  !GT?,P,!4<^ 

IMTFGFP  d  ,dmax  ,(j  .  v 

NFFD5  FO-T  in  T  n  RPCOVFR  'J  \  F  C  P  A  v  R  L  F  n  FOURIFP  CO^FF  I  C  I  FNT  c 

THIS  S'JRROUTINF  PFPLAO-F  x  j.  j  y  ov  jtf  FOURIFR  TRANSFORM.  ’-.'H EPF 
X(-)  +  IY(c)  -  TIJV,  t_o  ,  (  MOPTS-1  )  OP  (  X  (  T  )  +1 Y  (  T  ]  )  *Fy?  {  -P*T/ViCPTP  )  . 

PF&L  I  (  D  Y  A  X  )  y  P  (D'-OX]  y  C  \  D  *■  *  A  X  ,  D  ‘  ■'  A  X  }  ,  C  (DMAX,PMAX). 

•  a  (  (PMftX-n  **P  +  1  )  *  D  (  (  P'-'AX- U  **?  +  i  ) 

Rc  AL  T  (IF),  P  {  1  o  )  ,  r  (  p  ,  1  H  »  F  (  1-a  *  1  *>  )  *  A  (  1  'A  R  )  ,  n  MAS) 

PM  AX  =  ”1  F 

TV'0D  I  =6  •  7  F  F  1  P  S  3  0  7 
M=MOPTS 
lOn  cONTIN!t“ 

IF  (  Y  .  MP  .  (  /A  )  *A  )  r-n  Tn  A0  0 

PAfTOPS  °P  FOHP 

VO-  v 
V  -  V  /  A 

DO  FOO  J=1,M 

APG  =  T’FOD  t * p |_ n a  T  (  J-1  )  /FLOAT  (  vp  ) 
n=COF{  A  "  G  ) 

SI =SIN ( A  PG ) 
r  ?  =  ron  (  ?  .  n^-  a  DF  1 
S  P  =  S  I  A!  C  ?  ,0*AOG) 

C?  =  COF  (  F  .  a*ai?g  ) 
c^  =  PI<\|(q.n*Anr) 

DO  POO  <=MR  ,>'OcTR  ,‘*P 
J  1  =  J  4-  M  O 

J?=J1+M 

JF  =  JP  +  M 
J4=J3+v 

R1 =X ( J1 ) +X ( JF ) 

R?=X ( J1 ) -X ( JF ) 

1 1 = Y ( J 1  )+Y( J- 1 
T?  =  Y( J1  )-Y( JF  ) 

R  3  =  X ( J?  ) +X( JA  ) 

R A -X ( JF ) -X ( JA } 
n=Y(  JF  )  +Y(  JA) 

I a=Y ( JF ) -Y ( JA ) 

X ( J 1  }  =Fi+P3 

y  ( j  i  )  =  n  +iF 


non  non  on 


A. 57 


x(j?)  =  (R?4-iA)*n  +  (  r  ?-pa  )  *si 

Y  (  J?  )  =  (  I?-R4)*C1-(P?+IA)*-S] 

X  (  )  =  (  D  1  -P  "  )*C?+(  I  ]  -  r  -j  )  *  S  ? 

Y(J^)  =  (n-ii)*r?-(P'i_?i)#-3;> 

X(JA)  =  (o7-IA)*0‘5-*-(  T  7  +  o  tL  )  *  s  n 

Y  (  J4  )  =  (  I  ?+P  A  )  '3-  (  P  ?-  T  6.  )  *-$7 

?0o  c^TINUf 

300  CONTINUE 
GO  TO  100 
AQo  COMT  I VIJ- 

IF  (  v/jc  .  (  v/p  )  *?  )  00  TO  7  00 

FAfTOPS  OF  TWO 

MprV 

V=M/7 

DO  300  J=1 ,M 

ARG=TW0P  T*FLOAT  (  J-l  )  /FLOAT  ( vp  ) 
ri=ros(Ai?G) 

Sl=STN(ARG) 

DO  fCO  < = vp , mop T S ♦ vd 
JT  rJ  +  ^-^P 

J?=J1 +  M 

°  1  =  X  (  J 1  ) +X ( J7  ) 

P?  =  X ( J1 ) -X { J? ) 

I ! =  Y ( J 1 )+Y( J?) 

I ?=Y ( Ji ) -Y ( J? ) 

X { J1 ) =01 

Y  (  J  1  )  =  I  1 

X  (  J  7  )  =  p  7  *  r  i  +  T  7  *  S  1 

Y  (  J7  )  =  I  ?*C  1  -  R  ?  *  S  1 
300  CONTIMUF 

600  CONTIMUF 
GO  TO  A00 
70n  CONTIMUF 

IF  {  M  •  N  r  •  (  V  /  3  }  -3  )  CO  TO  1OO0 

factors  of  T h P H r 

MR  ~y 
M  =  V/F 

A  1  =  fOS  (  TWOD  I  /  -a  .0  ) 
oi=S  IN  (  TV,' DP  t  /f  ,  n  ) 

A?=COS ( 2  .  n*TV'0P  I  /p  ,  o  ) 

P7  =  f TM ( ? .0*TWOP  I /f  .0 ) 

DO  ROO  J=1 *m 

ARG=TV/0P  I  AFLOAT  (  J-l  J/-L0AT  (MR) 

ABSORB  TWIDDLF  FACTOP  INTO  ANALYSIS  COEFFICIENTS 

C?1  =  C0S  (  A-RG  ) 

S71 =S IN ( ARG ) 

C2?~  C?1*M-S?!*BI 


n  n  n  n  n 
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S??=C71*ni+S7’l#Al 

C??=C71*A7-S?l*n? 

S77=C71*°?+S7 1  *  A  7 
O1=C0S(?.0*ARG) 

SM=SIN(  2  •  0  *  a  R  G  ) 

007  =  0*  1  *  A  2  - s  -a  1  *Rp 
S  7  7  =  C 7 1  *  R  7  +  S  7  1  *47 
C 7 0 .  =  C 3  1  *■  A 1  -S'1!  *  R  1 
S‘5'»=CS1*^US'HA1 
DO  POO  <=vp  ,mopTS  »v,P 

J?=J1 +w 

JP=J?+V 

Rl=X ( J1 ) 

p  =  y( ji ) 

R7=x ( J7 ) 

I  7  =  Y  (  J  7  ) 

R7=X(JP 
I o  =  Y ( jo  ) 

X(J1)=R1+R7+Rp 
Y(J1  )  =  !1tI7+I'3 

X  (  J7  )  =  R  I  *  r  7  1  +  .  T  1  *  S  7  1  p  7  *>  7  +  T  ->  *  S  p  y  +  p  -J  *<-  7  -5  +  I  -a  *5;  -5  -5 

Y  (  j  7  )  =  t  i  *o  7 1 i  *  s?  i  4-  f  + 

X  (  J-3  )  =  R  ]  *0  7  1  -1-  f  l^S-n  +  3  1  n  I  7-^S^  J.'’  H  ^ 

Y(  jp  )  =  n  -17  1  *  S  -a  ’  4.  t  7*0p7-R7*SP7+  T  i«f  o'i-'J-JSS’'3 

POO  0  0  N  T I  V  M  F 
QO^  0  O M T I M I J  f 
GO  TO  700 
1^00  c  o  N  T  I  \! ! J  F 

IF  { M . M  c .  ( M/ p )  * P  )  on  TO  1  7  QO 
Faotorf  nc  vr 

VP  -  v 

V  =  v  /  7 

a  1 =oos ( t wo°  i  / 0 . 0 ) 

RT=S IN ( T WOP  I /p . 0 ) 

A  7  =  ^OS ( 7 ,o*TW0D [  /  p  .  n  ) 
o  7  =  S  T  >•'  (  7  .  'v*  T'-’QP  T  /e  #  n  } 

A  7  =  0  O  S  (  7  .  0*JU'0P  I  /  P  .  n  ) 

op  =  FIM(  7  .  o*TV'OP  I  /P  .  n  ) 

A  4  =  000  (  4  I  /P  .  "M 

P4  =  S  I  0!  (  4  «.  0  *  T V/ 0 P  ! /P .  0  ) 

DO  1 ?  0  0  J=l,v 

ARG=TV;0P  T -FLOAT  (  J-  1  )  /^LOAT  (MR  ) 

ABSORB  TWIDDLE  FACTOR  I  A! T 0  AMALYSIS  COEFFICIENTS 

C  7  1  =  COS ( ARC ) 

S?1=SIM( ADG) 
c  7  7  =  c  7  1  *  A  1  -  S  7  1  *  R 1 
S??  =  C71*R1  +  S?1'LAi 
C?^  =  C?1*A?-S?!*-E? 


S  ?  i  =  C  ?  \  *  R  ?  +  S  ?  1  *  A. 

C  7  Lx-  C  ?  1  *  A  ^  -  5  ?  1  *  R  "5 
$  ?  4  =  r  ?  i  *  n  ?  +  s  ”>  1  *  A  -5 
C  ?  c-  =  r  7  i  *  A  S  ?  1  4f  R  4 

S?  5  =r  ?  1  *  p  A  +  s  ?  I  *  A  4 
oi =r os ( ?.A^i°G ) 

Mi  =  STM  (  ?.o»APn  ) 

r  ?  ? = c  ■*  i  ■*  a  ?  -  s  ■=  i  *  r  ^ 

S??  =  C'3 1  -frRP  +  S'3  1  *■  A  ? 

5^^  =  n]^n4+.(;'!i  *.\l 
c  ?  Lx  =  r  ”3  1  *  A  1  -  C,  r?  1  *  n  1 
$  R  4  =  C  ?  1  *  °  ]  +•  S  “>  i  .  *  A  l 
r  O  R  =  f  7  1  *  \  n  _  C  7  1  *  p  ? 
^'jc-riixn^fC^uAi 
CA1  =CnS  (  G.n-^ARG  ) 

S  A 1  =  S  I  m  °  ,0*-APG  ) 
c  A  ?  =  C  A  1  *  A  ^  -  s  A  i  *  n  ^ 

R  A  ?  =  ^  A  1  *  m  O-  S  A  1  *  A 
r  6  -5  =  ^  A  1  *  A  1  -  R  A  1  v-a  ■) 

S  A  7  -  C  A  1  *  n  1  -*-  5  A  1  *  A  1 
C AA  =CA  1  SAA-.CAI  *RA 
Saa  =  tai  a  +  SA  1  *A A 
c  A  s  =  r  A  1  *  A  -5  —  S  A  1  *  ?.  ■? 

S  A  ^  =  C  A  1  *  n  ?  4-  S  A  1  *  A  7 

or  1  =cos  (  A  ,A-^AOG  ) 

SS 1  =  S  IN  (  A.O* ARG ) 
r5?=Cc'l“-AZi-Sc1"QA 
s  s  ? = r  *>  i  *  D  a + s  - 1  *  a  a 
*  A  °  —  S  c- 1 

$57=0^1  *  n  ^  +  3  G  ]  *  A  ^ 

C54  =  Ct?l*A2-Ss>  ]  * R ? 

S  R  A  =  T  ^  1  *  n  7  +  R  r'  1  *  A  ? 

C  *  c.  =  C  5  1'*  A  1  -  C  *  1  *  R  1 
c;  c ,  c;  _  (■  c,  i  *  o  -|  +  c  7  i  ^  A  1 

no  n  °r  <  =vp  »NOnT^ 

J1=J+V'-MP 

j?=ji  +»/ 

JA  =  J*a  +  ^ 

J5-JA+M 

R 1 =X ( J 1  ) 

U=Y(  Jl  ) 

P  ?  =  X  (  J  ?  } 

I ?  =  Y ( J  ? ) 

R^=X(J°  ) 

I ^  =  Y ( ) 

RA  =  X  C  J a  ) 

I  A  =  Y  (  J  A  } 

R'G^X  ( Jc  ) 

I  5  =  Y  (  J  0  ) 

X  (  J  1  )  =  R  1  4-0  p  +  9  o.  a.  P  A  4-  p  R 
Y  {  J  1  )  =  n+  I?4-!^4-TA4Tr- 


non  non  nnr» 
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•  P‘>*C?<?+I^*S?t5 

y  ( j*? )  =  r  i  *r?  i  -  p  i  *f? -i4.t?*-o??-3?*s2?  +  i^c?  ^-2^*0  ?^>  +  r  4  >0^-94*574  + 

xf^iroi-sc  01  +  11*0^1 +9?  *c',?  +  r?-*so?+9o*c33+i3*F'^+R/1.*rp4.  +  i4*cT>4.+ 

Y(  J?  )  =  J  1  *C  3  1  -P  1  *FP  1  -+-  I  o*C  ^  7-9?-*FA?-i-  [  G-*C  o  3  -Ro  *  ~  ^  3  +  I  4*0  tf.-94-k-.F'J  4  + 
.  TS*C^-O^*.50F 

X  (  JA  )  =  R  1  *C4  7  +  I  1  *0  4.  1  J-P  ->*r  4">  +  I  ?  *  34  7  J-~  'a*  C  4  1+  I  **  *  4.  9  A  *C4  4  4.  I  4*  3  4  a  + 

.  ,R^Ci^+  I  f-^Sff 

Y  (  J4  )  =  T  1  *CA  1  “9  1  *  34  1  +  T  ?  *C4?-9?  *  -4  ?4-  I  P  *C  4?  -RP*  340  4-  I  4  *  4  -*- 

•  if*c45-9c*-S4f 

x ( JF ) =Pi  -*r  f i  +  t  i  *sst+??*r^?+r?-^sc??  +  R'j*c:^  +  !'3*-r> 334-94*054+14  »■:' 5 4  + 
I  c  *3F  c 

Y(JF)  =  n*-^Fi-9l-*3F1  +  T?-*C5?-R  ?  *  35  ?->-  T  p  *f  5  i-Ri-K.^u[i(^c;A-:?4«5^  + 

•  J  c;  *rc;  e_o  e;  *5  c;  ^ 

i.l  no  contia1"- 
1  POO  CONTINUE 

GO  TO  looo  - 

1  30O  COVT  I  \U- 

I F  (M.Lr.  1  )  GO  TO  ?6nn 

0FNF9AL  factor f 

DO  1400  j=?vDvaX 
D=J 

IF  (  M.  po  .  CNV°  )  *D  )  GO  T0  1500 
14.00  COMT  I  vijc 

CALL  fct  frr 
i^n  C 0 N T  I  A ! '. .' r 

JT={P-1 ) **?+] 

cry  (JD  iPpjjDARY  FA:"T0P3 
DO  160  0  J  = 1  1 JT 

APG  =  TWpP ! * c L o  A T ( J-l ) / r L 0 A  T ( P ) 

A ( J)=COF ( A 90) 

P  (  J )  =  F  I  M  (  A  9  0  1 
16  00  COA’TINUc 

MO  —  *.* 

V  =  V‘/D 

DO  p 0 0 0  J='l,v 

ARG=TWOD I  AFLOAT ( J-l ) /FLOAT ( MR ) 

AP50RR  .T,-:iDnLr  FACTOR  IYTO  ANALYSIS  COFFFICIFMTS 

DO  1300  U=l»°  ’  * 

CHS  1  )  =  r  Of  (  FLOAT  (U-l  )  T-APG  ) 

5(U*  1  )  =  S  IN  (  cLOAT  (t!-1  )  *APo) 

DO  170  0  V=?  i  D 
JT= (U-l ) *  ( V-l ) +1 

C  (  U » V  )  =  C  '( U  » 1  )  *  A  '  J  T  )  -  3  ( U  >  1  )  #  o  (  J  T  ) 

■  S(U*VJ=C(U*1  )  *  ?:  f  J  T  )  +  F  n  J  *  1 )  *  A  (  J  T  ) 


1 


non 
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1700  CONTINUE 
lftOO  CONTINUE 

no  7  7  no  <=yr5»f!0DT^iMo 

G  n  M  F  R  A  L  A  M  A  L  Y  S  I S 

no  i oco  u=i*p  , 

JT  =  J  +  ^-'10+  {  U—  1  )  *V1 
R  (  U  ) =X ( J  T  ) 

T  (in  =Y  (  JT  ) 

|onn  ro^T  INt Jp 

'  no  7ioo  u  =  i [ , d 

X  T  =  n  •  ^ 

YT=°.n 

no  ? non  v  =  i  ,p 

XT  =  XT+R(V)»f(i!iV)+!  (VO-^M'I.V) 
YT=VT+  I  (  V  )  f '  »  V  )  -R  (  V  )  *c  ( '  U'M 
2  n  0  o  c  n  m  t  I  m  !  j  c 

jT=j+x-f-'D+  ni-i )  *v 
X ( JT ) =XT 
V( JT ) = Y  T 
7100  CONTINUE 
7700  COMTIMUF 
7 7 00  COOT  I  MU- 

GO  TO  1 ^ no 
7  4  00  COMTIMl.Jr 
RFTURN 
FM.n 


\ ' 


n  n  n  n  n  n  n  n  n  n  n  n  n  nnn  n  n  n  n  n  n  n  n  n  n  n  on 
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SUBROUTINE  Sopt  in 

UNSCRAMBLING  DPOGPiM  FOR  OME  0 1  v  p  N  S  1 0  N  A  L  COURIER  COcFFlClrNTS 

C 0 v-voN /CONST /NOP  TS 
COMVON/DAT/X  (  4P.0  )  ,Y  Ur  ) 

PFAL  S  (  4  S  0  ) 

IMTFGFO  JT 

INTEGER  00  *  L  I M  £  1  *=»  )  »  S T r °  (  1 S  )  »P  ,oVaX 

TNTfgfR  a ,p , r *  P  *  c  »F . g  ,M  » I i J »K ,L » AL  «  nL  «CL  *0L L  »FL  »GL *HL «  TL . JL ♦ 

.  <L  »Llr  *ML  *  AS  «  p  S  »  0  S  »  0  c  *  F  S  »  F  S  *  G  S  •  WS  »  I  c.  *  JS  •<?  *  LS  »WS 

r>  l  G  T  T  DrVFPScR  FOP  US^  VITu  POUR  in  .  5  ML'ST  Rr  THE  SA  ME  SI7C  AS 
x  a  m  n  y  . 

EQUIVALENCES  to  allow  indexing  to  s  ft  p  a  r  A  m  f  t  f  r  S  AND  ALLOW- 
scalars  fop  use  in  t h p  no  loops. 

POUT VALCNGF  ( as  ,  ptpp ( 1  )  )  *  {  PS  ♦  STfo  (  ;?  )  )  , ( rS ♦ STpd ( 3 )  )  , ( nS  *STFD ( 4 )  ;  , 

.  (FS,ST^D(S)  )  UFS,ST?d(a)  )  ,  (  g  p  *  S  T  p  p  f  7  )  )  ,  (  w.c  ,  S  T  pp  (  a  )  )  ,(  IS*STEP  )  }» 
.  (  JS.STfp  (10)),  (KS  ,ctpc  (11)),  (  L  S  •  S  T  p  D  (  1  7  )  )  »  (  N'  S  ,  S  T  p  P  (o)) 

FOU  I  VA  lc.NGF  (  AL  *1.  T  v  (  1  )  )  <1  (  nL  *  L  T  M  (  2  )  )  »  (  GL  *  L  T  V  (  7  )  )  »  (  ^  L  »  L  I'M  4  )  )  » 

•  (FLMT-  M)  )  *  (  PL  *L  I'M  A  )  )  *  (  GL  *  L  T‘-M  7  )  )  ♦  ( H  L  » L  I  */  (  a  )  )  ,  (  n.  ,  i  I  v  (  o  )  1  , 

•  (  JL  *  L  I V  (  l  0  )  )  ,  f  vr.L  »  L  T *v  (  i  l  M  »  (  L L.  <•  L  I  v  M  o  )  j  ,  {  ml  %  L  I  V  (  •  ^  M 

dm  AX  IS  S  rT  TO  A  o  P  F  p  WITH  pq!JR  in 
P.M A  X  =  1  p  ’ 

SFT  LIMITS  A  No  ST-p  ft/fs  ppom  j^ipp  LnoDS  GOING  OUT 
00=  1  7 

m= MOPTS 
100  COMTTNUF 

CHFCK  F0°  FACTOPS  0r  4 

IF  (  M.NF  .  (  M/4  )  *4  )  GO  TO  P 70 
M  =  V  /  4 

RFALLY  WANT  0-4*m-1  PUT  WE  GO  FROV  1  TO  4*M.  A.  fjppf  of  M  WITH 
MAXIMUM  DISPLACEMFNT  of  m  INITIALLY 

L I M ( DO ) =  4*V 
STEP ( DO) -M 
00=00-1 
GO  TO  100 
200  COMTINUP 

CHECK  FOP  REMAINING  FACTORS 
IF  (M.LE.l)  GO  To  son 
FACTORS  OF  2  »  F  ,  5  »  7  •»  11  » 1  B 


. 
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DO  ^00  JT=2 , ovaX  A. 63 

P  =  JT 

IF  (V  ,cf).  (m/d  )  *p  )  GO  jo  4QO 
700  COMTINUF 

C  • 

C  ERROR  r.X  I  T  Ic  FACTORS  4P,nVr  P'-'AX  N  F  c  D  F  D 

C 

CALL  FCT  ERR 
400  CONTINUE' 

MsM/D 

C 

C  REALLY  VM.MT  o  _  p  ^  v*  _  2  P!JT  v=;  ijpp  j  yo  d**<.  P  cjppc  qp  m  VI  TH 

C  MAXIMUM  INITIAL  D I EPl a CE vc\T  n:r  M 

c 

L I m ( DO ) 

S  T  F  p  (  D  0  )  =  *•' 
no=oo-i 
GO  TO  POO 
SOn  COMTINUF 

FINISH  OUT  THF  DO  LOOPS  TO  VA<F  OUTER  LOO^c,  FYFCUTF  ONLY  ONCE 

DO  G 00  JT  =  ]  ,00 
L  I  M  (  J  T  )  =  1 
STFD(JT)=1 
6 Op  CONTINUE 

SFT  JT  SO  THAT  JT  RIJNS  FROM  1  TO  N^PTS  IN  STFPS  OF  1  WHILE  M  WILL 
RUN  WITH  RFVFRSFn  DIGITS 


JT  = 

=  0 

DO 

7n0 

A  : 

=  1 

♦  AL 

,  AS 

DO 

700 

R: 

=  A 

_j 

a 

,RS 

DO 

jrn 

C  : 

:R 

,CL 

os 

DO 

700 

D  = 

=  C 

OL 

,DS 

DO 

700 

r  : 

=  D 

♦  EL 

*  FS 

DO 

700 

F: 

*  EL 

,FS 

DO 

7  00 

G  : 

,GL 

OS 

DO 

7  A  0 

H: 

=  G 

OL 

,HS 

DO 

7  on 

1  = 

=  H 

♦  TL 

*  I  S 

DO 

700 

J: 

=  I 

*  JL 

os 

DO 

7  00 

<: 

=  J 

OL 

OS 

DO 

700 

L- 

=  < 

OL 

OS 

DO 

700 

V: 

-L 

♦ ML 

♦  vs 

JT: 

=  JT+1 

S  (  J  T  )  =  X  ( ‘  M 
700  COMTINUF 

COPY  RACK  OUT  OF  THE  SCRATCH ’ ARRAY 

\ 

DO  BOO  JT  =  1  , NOP T S 
X  (  JT ) =S ( JT ) 

8 On  CONTINUE 


non 


C 
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JT: 

-n 

DO 

900 

A 

=  1 

*  AL 

♦  AS 

DO 

o  o  n 

0 

=  A 

♦  RL 

♦  °S 

DO 

9^0 

r 

=  R 

♦  CL 

♦  Co 

DO 

090 

D 

=c 

♦  DL 

♦  D  0 

DO 

ODD 

c 

-  0 

♦  FL 

♦  FS 

DO 

9n0 

F 

=  F 

♦  ~L 

♦  c  S 

DO 

9  00 

r. 

-  cr 

♦  CL 

♦  c-s 

DO 

0  90 

H 

=  G 

*HL 

♦  HS 

DO 

9  00 

I 

=  u 

♦  IL 

♦  IS 

DO 

900 

J 

—  T 

♦  JL 

♦  JS 

DO 

900 

=  J 

»’<L 

♦  <S 

DO 

900 

L 

-  / 

♦  LI. 

♦  LS 

DO 

990 

M 

=  L 

♦  *'L 

♦  '-AS 

JT  s 

=  JT  +  1 

S(  JT  )  =Y  ( *') 
ono  rOM  T  I  w,l  icr 

COPY  o A. r :<  OUT  OF  THF  SCRATCH  ARRAY 

DO  o^O  J  T  =  1  ♦  N  0  D  T  S 
v( JT ) =S ( JT ) 

9S^  CONTI r 
Return 
FND 


SUBROUTINE  FC  T  Rp 
FACTORING  ERROR 


A.  65 


C 

C 


C  FACTORING  FRROR  IN 

C 

C  CURRENTLY  T A < F N  IF 

C  ftRF  NOT  n I G  ENOUGH 

C 

V/RTTF  (6*100) 

FOUR  ID  OR  SORT  ID. 

A  FACTOR  A°OVF  13  IS  REQUIRED.  (THE  ARRAYS 

T^  mandlf  THINGS  ADQVr  13* ) 

100  -0RVAT  ( 1X*1 5 H FACTORING  ERROR ) 

RETURN 

END 

1 

\ 


FORTRAN  IV 


PROGRAM  TO  CALCULATE 
SYNTHETIC  SEISMOGRAMS 


n  n  n  n  n  non  n  n  n  n  n 
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C 


This  PROGRAM  CALCULATFS  SYMTHFSIZpd  GROUND  MOTION 

I 


120  FORMAT  (  2 F6 . 1  »  2.  I  G  ) 

121  FORMAT  (  1  HI  »  2  HlJN  !  T-  ^6 . 1  » IX  ,  AHNDT  — *  !  6  » 1  X  ,  3H  I  P  =  *  I  6  ) 

122  FORMAT  (  1  X  «  2FHF  NO  OF  THF  RUN  »  T  H  A  N  <  YOU) 

122  FORMAT ( 1  HI ♦?5HFRFOUFNrY  OF  THF  PULSF  IS  * ^ . 1  *  1 X » 2HCPS ) 

TP...NUIMRFP  OF  FFfOMnS  PfD  INCH  TO  PL°T 
Pf'6..,du(S-  FRFOUFN'cy 

COMMON /  VFR T  /  pWr)TP  (  'tn^n  )  /hOR/PUdTP  (  ?onn  )  /T  jM/pcR  (  ?onn  )  /  fR/fRfQ.(  300n 

.  )  /  R  ^  A  L  /  H  R  P  (  0  O  C  )  *  7  R  p  (  ^  P  n  0  )  /  1  U  a  F  /  H  T  P  (  F  n  2  )  ,  7  |  D  (  Q  o  n  2  ) 

COMMON /  SC  A  L  F  /  H  H  A  X  , VM  a  X  /OAT  /X  (  aAAo  )  ,Y  (  ^on  )  /CONST /NOP  TF  . 
COMMOM/CST/ANG1P  »UN  I  T  ■  PFR  *N0F  ,.NPT  9  IP  jCINC/ARfa/pUF  (  F-on  ) 


RCAn ( 5  9 1 ?0 )  UN  IT  »PCR  9  NOT  , IP 
W  R  I T  F ( G  » 1 2 1  )  UNIT,N°T,IP 
WRTTF(6*1?2)  pFR 
NOPTS=NpT 


K*r=4*S00 

CALL  PLOTS ( PUF (  1  )  »xr  ) 
CALL  PLOT(A.0,?4.n,-i  ] 
c 

L  =  1 

CALL.  TRANS  (  L  ) 

C 

A  I  =  (  FLOAT  (  N°T  )  /UN  T  T  )  4- 1  o  . 
CALL  P L O-T  (  A  I  tO.O  .-"5  ) 

CALL  PLOTfo. n,o.o ,ooo ) 

WR I T  F ( 6 i 1 2  2  ) 

STOP 

CND 


. 


, 
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SUBROUTINE  LAGS(m,C) 
ppal  C(!n) 

KsrV/7  +  l 

L=M/?-l 

no  i  1=1 ♦  !_ 

C( i+k )=r ( <- i ) 

1  rONTIN'LJc 
V/R  T  T  F  (  ft  *  ?  ) 

7  FORMAT ( IX  * lOHAFT^p  LA-S/) 
RETURN 
rND 
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SUOROUTINF  rORR  (MO‘Sc  I  ) 

RF  A  i_  mqm»FI 
I  P ( )  1  » 1  » A 

1  T  F  (  F  I  )  ? 

?  A  M  G  =  0 . 0 
GO  TO  7 
0  A  N  G  =  1  •  0 
GO  TO  7 

A  T  F  (  p  T  )  F  »  *  *  A 
*  A  A]  G  =  —  1  •  0 
GO  TO  1 

6  A  M  G  =  0  #  0 

7  mom  =  MOM+  ANG*'5  •  ]  4  1  F  °  ?  F  K 
RF  TURN 

END 


SUOROUT I  ME  MAX ( Mix »Xv AX  ) 
DIMENSION  X ( 1 0 ) 

XMAX=0. 0 
00  ] 0  1=1  » M 

IF ( APS ( X  (  I  )  )-XMAX )  10*10*11 

1 1  xvax  =  ars (x ( t  n 
10  CONTINUF 
WR I T  F ( 6  * ? o ) 

? 0  FORMAT (  1 X  *QHAFTEp  max/) 
RETURN 
FMO 


SUBROUTINE  FAC  TOO(X»N') 

Rfal  X ( 1 0 ) 

no  10  ! = 1 ♦ N 

X (  I  )  =  X (  I  )  /FLOAT { w ) 

10  CONTINUE 
RETURN 
FND 


SUBROUTINE  PPR  IOr>  (  p  ,m  ,  TF  ) 

RFA|_  P(10) 

DO  ]0  I = 1 *N 

P (  !  } = FLO AT (  I -1  )  / ( t REFLOAT ( N )  ) 
10  CONTIMUP 
RF  TURN 
FND 


■ 


SURROUT  I  NF  COMCON ( V » N ) 
REAL  V ( 1 0 ) 

DO  10  I = 1 
V (  I  ) =  -V (  I ) 

10  COMTINIJF 
RFTtJRN 
FNO 
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SUPROUT  I  NF  PULSE  1 

COMMON/ VF^T  /Dv/DTP{FOnP)/HOR/PMr}TP(Fnpn)/TIM/DPR(?POP)/FR/FRPQ('2PPn 
•  )/RrAL/wPP(^non)  ,  7  p  °  (  1  o  p  n  )  /  1  N  A  0  /  H  T  D  .  7  T  D  f  ^  p  n  o  ) 

COMMON  /  SC  ALr  /HM  AX  ,VMAX/OAT/X  (  )  ,y  (  )  /  cON  c-  T /NO  P  T  F 

COMMON /CRT/  AN  01  p  LIN  IT  »  P  FR  »  NO  r  « N  D  T  >  !  p  » 0  I  M  C  /  A  R  c  A  /  P  ’ F  i  5  0  0  ) 

DO  1  1  =  1  » NOP 

IF  (PFP-FRF0(  I )  )  7  ,7 

*  A  =  l. 

GO  TO  4 
7.  4  =  0. 

4  HPP ( I ) =NQD ( I ) *A 
HID(I)=HIP(T)*A 
7  P  P  (  I  )  =  Z  R  P  (  I  )  *  A 
Z  I P  (  I)  =  7  T D  (  I  )  *  A 
1  CONTINUE 
VRITE(6*5) 

F  F n p w  at  (  1  X  » 1  ? H A F T <" p  PULSF1  /  ) 

P  F  T 1. 1  p  N  i 
PNO 
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SU° ROUT  I Nc  TPLOT(Til) 

rOWON/VFRT  /  Dv/nTD(?^r>n)/HO^/P:inTP(^onn)/TiM/PcrR(^nno)/FR/(:rRF'Q(?onr1 
,  1  /PrftL/'-'RP  (  )  ,7RP  (  ^non]  /[Vag/mtp  (3n.n0)  ,7  jp(  ) 

COMMON /SC A LF/HMAX , VM A  X/DAT/X ( S OOO )  ,  Y ( 3 POO ) / COM S T /NOP T S 
COMMON /CST/ANO  IP  *UM  T  T  »  opp  *NOF  ,NPT  » IP  *G I  NC/AP.F  A/PtJF  (  5  on  ) 

CH  =  -1  .  S 

X  S  C  =  1  .  /  (  T  P*IJN  I  T  ) 

HSC=? . /UMAX 

v  s  o  =  ? » /  v  m  a  x 

MD  =  fU  D  J  /  2 
TUM  I  T  =  UN  1  T 
DO  i  1  =  1  » N  p 

CALL  PLOT  (  FLOAT  (  I- 1  )  #  X  S  C  * P W D T D  (  I  )*VSC  *2  5 
1  C  O N>  T  I  m  IJ  F 

CALL  SYMP.OL(P.l  »?.?»o.lS  j^PHVcRTICAL  GROUND  MOT  I  ON  »  0 . 0 , 2  2  ) 

CALL  DLO^ ( .0 , -A. 0 ,_7 ) 

DO  ?  1  =  1  *N,P 

CALL  PLOT  (FLOAT  (  1-1  )  #  X  SO  » D  '•  J n  T  D  (  I) *  u  S  C  •>  ?  ) 

?  CONTI  MU'7 

CALL  SYMBOL  (  P  .  1  »  2  .  2  <*  0 . 15  »  ? 4 H H 0 R  T  Z 0 N T  A L  GROUND  MO T  I  ON  ,  n  .  0  ,  ?4  ) 

GALL  DLOT  (  n  ,  0  ,  CH  ,-p  ) 

CALL  PLOT  (n.'i^CHi-1  ) 

CALL  pympQL  (?.nO,n.ntiD;,i  ?HORi  !I^taL  ‘-'ODFL  .n  .  n  ,  \  1  ) 

CALL  SYMROL(2.0’-0.7* ".1 S  »°SHANGLF  OF  INCIDENT  DEGRFFf, 

.  0 . 0  ,  ?  P  } 

V!o  T  T  p  (  A  ,  A  ) 

C  A  L  L  N  U  M  p  P  R  (  A  .  6  *  -  0 . 7 , 0 . 1  F  «  A  N  G 1  P  ,  0  .  n  ,  _  i  ) 

CALL  SYMBOL  (?.  0  *~1  .?*  15  »  ?  1  HPULSP  FR^Ol-PNCY  IS  C.P.S.* 

.  P  .  ^  ,  ?  1  ) 

WR I T  E < G  » 4 ) 

CALL  WIJMPER  (  A  .  A  , -1  > 

CALL  dLOt ( o . 0 , -6 .n ,-o ) 

WR I T  F ( S  *  A ) 

WR I TF ( A , h ) 

u  format  (  1  X  » 1  p HTPLOT  I  c-  OKAY) 

RETURN 

fmd 


. 
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SUBROUTINE  TRANS (L) 

C . . . . . . . . . c . 

c 

C  CALCULATION  0  c  T  R  4  N  c  r  f  r  c ! <  N  f  T  I P  v 

C  APPARENT  surface  VELOCITY  VUSt  RE  GREATFP  THAN  P  VF  LOC I  Tv  IN  ThF  TOP 
C  LAYFP  A NO  S  VFLOCITY  IN  A N v  LAYER 

C  i 

•J  OF  Foomat  (  y  T  6  ) 

1  0  1  FORMAT  (  t  ,tr1  ‘i,4  ) 

1C?  FORMAT  (  1XmOHTHIC<./yy/»?X  >  1  4HP  VE  L  . /KM/SFC/ - 1 X  *  14HS  VEL . /KM/SFC/ • 
.  ?X  > 1  >.HOFv  .  /r.o/ff'r/  /  ) 

1  0*>  F  0  R ' '  A  T  (  PF1  0  .4  ,  IE  ) 

1  04  FORMAT  (  ?X  »  OHAD  .  V-  1.0  r  .  ,  4X  9  :  .  TNC.aOT.  >  4  X  ,  1  ^  u  A  N  G  •  INC  .  SURF ,  ) 

105  FORMAT  {  r  1  0 , 1  »  ,  7  X  .  r  1  ?  ,  *>  ,  1  X  ,  c  1  ?.?*‘aX«cl?.‘a) 

106  FORMAT  (  ].X  »5HLAYFR  ,  RX  ,  1  1HTH  rr<  RATIO  ,?X.,7HP  RATIO, ?X,7HS  RATIO, 3X, 

•  1  0  H  0  F  v  .  o  a  t  I  o  /  ) 

•107  F0RMAT(1X,T?*6X»“!2.A*‘3^10.4/1 

I  OR  FORMAT  {  1  X  *F7  .  1  ?F7  .  fv ,  a  R7  .  ?  ,  r<  r  .  tv ,  ->X  ,-R  .  4  ,  7X  ,  ?F  6  •  4  ) 

1 0o  FORMAT ( 1 X • 1 ?HcN0  OF  TRAMS) 

no  FORMAT  (  1  HI  ,  1  7H  TRANSFER  f;  Jnc  T  I  ON  /  /  ) 

111  FORMAT  (  1  X  »  6  H  D  -  P.  I  C  0  i4X  »  4  H  r  R  e  Q ,  4  X  >  2  M  T  W  >  5  X  ,  ?uopy  ,4X  ,?HTU  »  5X  ,  *a  HP  HU  »4X  , 
.  4 7 HHOR  T  M ZOM  T  A  L / R R  A  L - I m  A  G / , VF  R  T I C  a L / R -  A L  -  T'  A  G / /  ) 

117  FORMAT  (  -*x  ,  4  HANG  b  7X  ,  "  HO  A  NO!  >  5  H  N  n  A  MG  /  ) 

i  i  ^  FORMAT  (  X  ,  4HG m  I  N  >  GX  «  4 U'G  I  VC  ,  7X  .  r»HMOr/  ) 

II  4  c  O  R  M  A  T  (  i  X  *  4  H  /  /  ) 

1  1  6  FORMAT  (  1  HI  ,  nHCRUSTAL.  mqofl»  !*-//) 


rOMMQM/V-PT  /  P  W  0  T  D  (  *>000  )  /uOR  /  °!  ,r>  T  o  (  iAAn  )  /  T  TM/pco  (  )  /  F R  / F  !?FQ  (  ?  070 

,  1  /  R  f  A  L  /  H  p  °  (  *>  o  o  o  )  .  7  R  P  (  T'^O)  /  I  M  A  G  /  h*  I  p  n  n  ~  0  )  •  .?  T  p  f  3  nn  ^  ) 


COMMON  /  SC  a  L  -  /mm  A  X  «  VY  AX/OAT/X  (  "GAa  )  « Y  (  7^0^)  /  "ON  ST  /NO0  IF 
COMMON  /  CS  T  /  A  *1 G  1  D  ,  l  ;»■'  I  T  »  o  FO  ,  Mn  r  ,  _v  DT  ,  T  O  ,  r.  I  VC 
n I  Me NS T  ON  0(60)  » A  (  ~  o )  ,  P  ( 6 ^ )  ?RmO ( 5° )  ,  OGAM ( c R  )  , AV^  ( =0 ) 
. ORP ( 50  )  , OH ( 6  0  ) 


A.  (60) 


1 


CALCULATION  O^  TRANSFER  rUNC  T I  ON 
RFAO(6»l00)  NC'M 
DO  7  6  I  A  I  =  1  «  NOM 

9PA0  T  H  c  T  Oc-'NT  T  F  T  C  a  T  I rNV'  OF  Thc  CRUSTAL  MOOFL 
wOl  .  .  .MljYFFp  nr  LaYFDS 
P  F  A  0 ( 6 , 1 0  n )  MOL 

rfao  THF  LAYFRS  PARAMETERS.  o  =  thi  cknfss  ,  A.=p  VELOCITY,  o  =  S  velocity 

RHO  =  Oc:n r  r  t  Y 


WP  I  TE ( 6  >  n  6 )  I  A  I 

WRITE (6 >10?) 

no  R  I=1*M0L 

R  f  a  0  (  6  »  1  0  1  )  0(  I  )  sA  (  I  ).^D  (  I  ]  ’  RHO  (  I  ) 
WR I T  F ( 6  >  1 0 1  )  0  (  I  j  »  A  (  I  )  , p  (  I  )  »  RmO (  T  ) 
?  CONTINUE 
WRITE (6 >114) 
wr  i  t  f  ( 4  >  n  4  ) 


non 
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M  =  mol 

Ml =MOL-1 
N?=MOL-? 

C 

C  R  F  A  n>  AMG1  =  I  N  I  T  I  AL  A  N  G  l  c  or  IMOIO^MOF  at  TH^  SWRfacF 

c  r'AKjrjrTMfRFMFMT  op  J  M  c  AMOLrS  OP  JMf  Tr>rWfp 

C  MOAMG^NU'-’PFR  OF  ANGLES  of  INCIOFNC-  TO  °F  OALv"UL  at  - O 

c 

RFAO( 5  » 1  OF )  A  MG  1  »  D  a M G 1  ,  NOA.MG 
WR  I  T  F  (  6  »  11  ?  ) 

WR  I T  E ( G  9 1 0  7  )  A  MG  1  , D  A  M  G 1  *  MO  A M G 
WP I T F ( G  9 1 1  A ) 

WR I T  F ( 6  *  11  A ) 

c 

c  Gvi  I  M  IS'UV  FRQIJFMCY  O^GI^FO 

G  I  Mf  =  I  NOPF^RNT  or  FPfOUFMC  Y 
MOF  =  MUv'FFR  of  VAL.U-S  Of  r  p  r  q  ij  f m C Y  OFSIRpo 

RFA0(F,10^)  G^IN»GINr »N0C 
WPTTF(G»11 

WR  I  T  F  (  p  ,  l  ^  7  i  G  V  I  M  9  G  I N  0  »  NO' r 
WP  I  T  F ( G  9 1 1 4 ) 

WP I T  F ( 6  » 11  A ) 

•  A  M  G 1  P  =  A  m  G 1  ~  0  A  M  G 1 
00  A  T  A  =  1 »  N O A  M G 
A  M G 1 P  =  AMGi P +  0 A MG  1 
S  I  M  I  =  G  I  M  (  A  M  G 1.  P  /  F  7  .  ?  o  f  ?  q  ) 

S  T  M  T  m  =  A  (  1  )  a-  p  I  A'  I  /  a  {  MOl  ) 

0- A ( MOL )  /  S  I  M  I 

A  I  V  -  A  D  s  I  M  (  S  I  A'  I  Ml  *  F  7 . 9  O  C.  7  Q 
WP T  T  F ( G  i 1 O  a ) 

WR  I  T  F  (  G  •»  1  0  F  )  0  »  A  MG  1  P  «  A  r  W 

WR  I  T  F ( G  9 11  A ) 

WR I T  F ( 6  *  1 1  A ) 

00  f  m=L,MOL 
ro\/A  =  r/A(M) 
rovn=r/o ( m 1 

OG  A  v  {  w  )=?./(  o  OV-  *«-  7  ) 

OGAO/i  (  v  )  =06 A v1  (  w  )  -i  . 

ORA ( ^ ) =G°PT ( ARS ( roVA*#?-l . ) ) 

OPP  (  M  )  =  SORT  (  APS  (  OOVP-&  -y-p-l  .  )  ) 

OH ( V ) =  P  H  O ( w ) *6# C 
f  GOMTINUf 
TH=r.n 
00  G  L  =  1  9  M  ] 

TH=TH+0(L) 

G  GOMTIMUF 

WR  I  T  F  (  G  9  i  p G  ) 

00  7  1  =  1  » MOL 
PK  =0 (1  ) /TH 
OK  —  A (  I  ) / A (MOl  ) 

VK  =p.  (  I  )  /  P  (  MOL  1 
RK=RHO ( I L/RHO ( MOL ) 

WR  I  T  F  (  G  9  I  0  7  )  I  9  P  '<  ,  Qv  9  V’<  •  R  K 
7  COMTIMIJF 


- 


WR  I  T  F ( 6  *  1 1 0 ) 

WR I T  F ( ft , 1 1 1  ) 

no  p  I F p  =  1  » MOF 

G  A F  =  F L O  A  T (  T  FR-1  ) *G I  NT 
FPFO (  ! FQ  )  =  G A F 

!F(rDtrO(  IFD)  '"■n  Tn  \ 

opp( 

]  WVN0=6. RFM pFF*FRcq(  i  fr  )  /  C 

A]  1  =1  . 

A  1  R  =  0 . 0 
A  ?  1  =  0  # 

A  R  ?  =  1  . 

A  P  1  =  0  . 

A  p  R  =  0  . 

A  4  1=0. 

M?rA. 

C  0  V  P  U  T  r  F  L  F  w c  m T  S  O  R  A  MATRIX  ^OR  LAYUPS 
no  o  n*  =  i  *Nl 
g  a  n- nr-  a  m  ;  v ) 

GAMM1 =OGAmi (m ) 

R  A  =  0  R  A  (  V  ) 
o n  =  P) p  n  (  k*  ) 

H= HM(V) 

P=mVMO*D  ( )  *  R  A 
O=WVN0*D  ( K‘ )  x-RR 
S  I  D  =  F  I  N  (  p  ) 

W=5 I  P/p  A 
P  P  =  P  A  *  F  I  P 
rOP=GOG ( D ) 

S  1 0  =  S  I N  (  0  ) 
rp=fio/pr 
Z  =  P  p.  *  G  I  0 
roo=cos ( o) 

PHO—PHR  c  ) 

n  v  =  n  (  M  ) 

R  1  1  -  G  A  m  *  C  0  °  -  G  A ' 4  ‘A  1  -x-  f  O  r' 
p  1  ?  =  G  A  n  1  x-'y  4-  G  A '  ■  "■  Z 
Dp  =  _(ror-r;AO)  /H 

P  1  A- ( W  +  Z ) / H 

071  =  +  (  G  A  D  +  G  A  M  1  *■  R  p  ) 

P  R  R  =  -G  A  M  A1 1  x-r  O  p  +  G  A  'A  *  c  P  n 

07^  =  _-(pp  +  opi/f-i 

p  R  4  =  R  1  P 

p  7  i  - f|x-G  a  x- g a M'i  i  x-  (  GOP-r 00  ) 
p  p  7  _  u  x-  (  G  A  ‘A  .V  1  *  G  A  V  */  1  x-  V/  +  G  A  A’  *  G  A,'-'  *  Z  ) 

R  P  p,  -  p  7  7 
P,  P  A  =  P.  1  R 

R41  --H#  (  GAY*GA‘A*PD+G a  V'M  XGA‘A'11  *RR  ) 

R  A?  =  °*1] 

R  A  p  -  R  p  1 
04 4 -PI  1 

r  A  1  1  =  P  1  1  X  A  1  1  +  o  1  7  X-  A  R  1  +  P  1  p  *  A  P  1  +  R  1  4  x-  A  4  1 
F  A  1  R  =  R 1  1  *  A  1  R  +  °  1  ?  -y-  A  p  r  +  n  ]  p  »A°  ’  J-c>  1  4  *  A  4  R 
PARI  =  n  ?  ]  y'  A  1  1+QRRXrARl  4-  R  p  P  *  /  p  1  +PRAXA4  1 


' 
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F  A  7  7  =  R  ?  1  X  A  ]  7  4-  P  7  7  X  A  ?  7  4-  p  7  0  *  A  P  7  4-  p  7  4  *  A  4  7 
F  A  7  ]  =  P  P  I  *  A  1  1  4-  R  7  7  *  A  7  1  +  P  p  °  *  A  t  1  +  o  -3  4  *  a  4  1 
F  A  7  7  =  R.  7  1  *  A  1  7  4-  P  p  ?  *  A  7  7  4-  P  7  Q  *  A  P  7  4-  d  p  4  *  A  4  7 

FA41s=P41*An4-P4?*A?]4-P4'a*AP'L+o44*A41 

P  A  4  7  =  p  4  1  *  A  1  7  4-  p  4  ?  X  A  7  7  4-  P  4  p  *  A  7  ?  -u  o  4  4  *  a  4  7 
A  1  1  =  F  A  1  1 
A  1  7  =  F  A  ]  7 
A?1 =  F  A  7  1 
A  7  ?  -  F  A  7  7 
API =EA?1 
A  7  7  =  r  a  p  7 
A  4  1  =  P  A  4  1 
A  4  7  =  P  A  4  7 
o  roNT  im;ip 
•  A  7  1  -  —  A  7  1 
A  4  1  =  -  A  4  1 
Gi^rDGA’^  { N  ) 

GAM''1 1  =PGA  v,]  (  \i  ) 

R  A  =  0  P  A  (  N  ) 

P  p  =  p  P  R.  (  A]  ) 

H=  OH  (  Kl  ) 

C.  GOMDtJTP  FLFMFiNTS  of  p  P'VP^S-  FOP  THr  LAST  LAYPP 
p  ]  t  =  _  r3  a  1  ‘  x  c  P  V  A  *  7 
P 1  p  =  1  .  /  (  R  HO  (  N  )  *  A  (  M  )  -x-  A  (  v  )  ) 
p  7  7  =G  amm  i  x c 0 V A  -xx  7  /  R  a 
P74=p IP/PA 
P44= 1 . / ( MXG AM } 
dpp  =-Q44 /P° 

Q  -3  1  =  _  P  P  P  X  G  A  1 1  ‘4  1  x  H 
p  4  7  =  I  •  0 

P  A  n  =D  1  1  *A.  1  1  4  P  1  0  x  A  P  1 
P  A  1  7  =  P  1  1  X  A  1  7  4-  P  1  P  X  A  P  7 
F  a  7  1  =  P  7  7  x  A  7  1  4-  P  7  4  -X  A  4  1 
F  A  7  7  =  P  9  7  *  A  7  7  4-  P  7  4  x  A  4  7 
PAP1=PP1#A]  14-DpPXA.P1 
PAPP  =  PP1*A]  7  j-  p  7  7  X  A.  p  7 
P  A  4  1  =  P  4  7  X  A  7  1  4-  P  4  4  X  A  4  | 

F  A  4  7  =  P  4  7  X  A.  7  7  4-  P  4  4  X  A  4  7 

P  R  =  F  A  7  1  X  P  A  77-PAT  1 x  P A  4  7 - P  A  1  p X  P A 4 1 -t- f A  7  7  x  p  A  p ] 

P  T -PA  1  1 XP A  7  74P A  7 1 *P A  4  7 -PA  1  7 XP  A  7 1 _C  A  7  7XP A 4 1 

D  F  N  S  0  =  U  R  *  D  P  4-  R  I  xp  \ 

UP  M R  =  E  A  7  7 *-R  j  _ p  a  4  7  *0  P 
UPN  I  =PA  d  7xnP4-r  A  A  7  XR  I 
WPN I =F  A  4 1 xpo  +  P A  p ] X n I 
V./D  \t p  =  _ F  A  P  1  X P P  +  F  A  4  1  X 0  I 

PV/DTP  (  T  F  ?)  =  ((  7  .  /  PPMPO  )  xpQP.T  (  VD  N  R xm  p m R 4-',/PM  I  #yPN  I  )  )  *C OVA 
PPHWP  =  ATAN'  (  +V/PN  T  /WPNR  ) 

PUDTP  (  I  p R  )  =  (  (  7  .  /HP  MSQ  )  *-5.QR  T  ( I  !DNR*  1 IPKR4-UPN  I  *IJP«NT  )  )  #C0VA 
PPHUPsA  Tam  (  4-UPM  I  /UP  Y  0  ) 


■ 


n  n 
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C 

r  A  LI.  ro^R  (  PDH’/,'d  *  l.\,PM  T  ) 

TALL  CODD  (  °PHUn  *U°N|  I  ) 

HI  P  (  IFR  )  =PlinTP  (  IF!?  }*MN  (PPHUP  ) 

HPD(  IPR)-PURTP{  I  FP  )  #FOS  (  ODHUP  ) 

Zdd  (  I  F R  )  =  =>-WDTD  (  I  c  p  )  *r^P.  (  DPHWP  ) 

Z  I  D  (  I  FD  )  -D'.fP)TD  (  I  -F  P  )  -M-  ^  I  .^J  (  PPHVP  ) 

c 

«*:R  TTP  (  A  >  1  OP  )  P  CR  (  T  "R  )  ,fRcC!  T^R  )  » Pv.'RTR  {  \  rR  )  ,ppw,;d  ,ci|nTO(  |FR)  » 
. PPHUP »HRP (  I  PR )  » H I D (  I FR  )  ,  ZRP (  I  c R  )  ,  Z  T  D ( I c R ) 
o  CONTINUE 
c  . 

IF(L.FO.l)  call  mOTION(IAI) 

IR(L.cO.P)  C  ALL  TDI  !Mf"  (  T  A  I  ) 

c 

4  r  n  m  T  i  m  u  f 

?  F  CONT  I  RU  p 

WR T  T  ^ ( 6  » TOO) 

RETURN 

FMD 


* 


no  n  n  n  n 
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SUOROtJTINF  MOTTON(IAI) 

c . 

c 

COMMON/  VFRT  /P''fnTP  (  ^  nrn  )  /hCR  /  Pi  in  TP  (  ?  oro  )  /T  I  iM/DcR  (  Rooo  )  /FR/FRoQ  (  soon 

COMMON  /SCALP/  u'-VX  ,VvaX/D4T/X  (  )  .Y  (  A  A0A  )  /CONST  /  N  0  P  T  S 

COMMON /CST/ AN G 1 P  *UN  IT ,OFR , NOF ♦NPT • TP , C I MC / a P F * / RUC ( 500 ) 

c 

117  FORMAT ( 1 X,4( T5»F]?.4*1X»C1’,^) ) 

HR  FORMAT  (  1  HI  »AX  ,  HMT  I^F-V  motion//) 

1 T O  FORMAT (  UH  » SX H 7HT T M--M  MOTION//) 

170  FORMAT!  1  H]  »  RHHMAX-  •  Ff  #  4 , 1  X  i  Fh'VM  AX-  «  CG . A  ) 


CALL  DULRF1 


N  i  =no°TS  /? 

CALL  7.FR0  {  NODTS  ,  X  ) 

CALL.  ?FRO  {  NOP.tf  ,  v  ) 

CALL  VC'/F  (NOF  ,7RP  ,  X  ) 

CALL  mqvf f  NOF , 7  I p , v ) 

CALL  LAC-f  (  nop  tf  ,  y  ) 

CALL  L  a  o  S ( M  n n  T  F , V 1 
CALL  COM CON ( Y  »  N 1 ) 

CALL  FOUR  10 
CALL  Sooj  TO 
CALL  FACTOO ( X *N0DTR ) 

CALL  FACTOO  (  Y  ,  NODTF  ) 

CALL  PFR 100 { nFR »NODj f mo ) 

C  A  |  L  MOVF  (  MOOTS  *  X  *  D!OTP  ) 

CALL  Z  f  P  n ( M  0  D  T  S  »  Y ) 

FALL.  7ron(MrPT.c,Y) 

CALL  m  o  V  F ( NOF ♦ H  o  D  » X ) 

CALL  MO V F (NOF iH!niY) 

CALL  LAGS (NQDTS ^X) 

CALL  LAGS (MOOTS ,Y) 

CALL  C0MC0M(Y,M1) 

CALL  cOi  )R  10 
call  SORT  10 
CALL  FACTOO ( Y, MOOT F) 

CALL  FAHOO(XiNOotf) 

CALL  mqvF (N0PTS »X ,DUOTo ) 

WR  ITE ( 6 » T 1 9 ) 

WR  I T  F ( 6  *  1  17)  (  I  *PcR(  I  )  »O|i0TD  (  I  )  »  T  =  1  iNpT  ) 
WR  I  T  F ( A  ,  n P ) 

WR  I  T  f  (  5  ,  n  7  )  (  I  *  D  r  R  (  !  )  9  o'<;oTD  (  n  ,  1  =  1  ,Ndt  ) 


r>  n 
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CALL  VAX (NPT *PUDTD ,H^AX) • 
CALL  MAX  {  NP  T  »  P  W  DTP  »  V1-'  AX  ) 
V\'  R  I  T  F  (  6  »  1  2  0  )  H  M  A  X  *  V  A  X 


CALL  TDLOT {  I  A  I  ) 

RETURN 

ENn 


FORTRAN  IV 


TRUNCATED  TRANSFER  FUNCTION 


PROGRAM 


n  n  n  no  n.n  nnnnn  non 
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this  program  calculates  tpuncateo'  transfer  function 


oo  FORMAT  (  l.X  »  P^HFNP  o  F  tuf  pum.  thany  yon) 

10  Ft  cr)PM/>T(?lft,^p7,pj 

10 1  format  (  1  H]  *6HN0?TS-  •  I  6  *  TO-  *  T6  •  1  X  OHCFR  .Ffi  .  i  ,  ?  ht-  •  r  6  .  ?  *  *  Hi  IN  I  T-  • 

•  F6 .2  ) 


I G  •  •  •  NT  JMQ  PR  0C  SDECTRA  POINTS  TO  PLOT 
MOP  TF  .  .  .  I N  T^O.p  A  T  I  ON  L  T  M  T  T 
T.  .  .TRUNCATION  TIME 

COMMON  /  VFRT  /  PHDTP  (  i  n  p  ^  )  /HOP  /  P>.  'HT  p  (  -a OP  A  )  /  A  Rf  a  /  oup  (  F  p  ^  5 
rO^MON  /OAT  /  X  (  -aopp  )  «Y(  -a^OO)  /  T  i  m/PfR  (  )  /rR/p^po  /  )  /rON  ST  /NOP 

rO^MON/PcAL/HRD  (  ?00n  )  ♦7RO(^r»OOJ/TMAr1/HTP(70or'J  ,7  fP  ( 
COMMON/SCALF/H  M  A  X  ,  V  * 1 A  X  /  c  I  N  r  /  0  r ]  c 
CONWQN/ROL/ A.NG1  P  *IJN  I  T  *  IG  »NOF  ,f  FR  *  T 


RFAO  (  5  » 1  OP  )  NOPTc  .-TG  .  <^fR  ,  T  .1  IN  T  t 
'•fp  I  TF  (  £  ,  1  n }  )  M0Dt  f  ,  Tf-,rcRiT  Y'MIT 


L  =  p 

CAl.L  TRANS  (L) 


CALL  PLOT(  10.0*?."  «-**  ) 
CALL  PLOT ( o  .  o ,p . ^  ,  poo  I 
V'P  I  TF  (  A  »  °°  ) 

^T  0° 

END 


* 
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GURROUTTNF  T  I  *-4E  (  I  A  [  ,  <  .  L  ) 

C.OMMOM/Vr  RT  /PUD  T °  (  oopo  )  /MOR/Pi.'DTP  (  poop  )  /aPfa /3I  'F  f  5^n  } 

C  O^'MON  /  R  F  A  L  /  H  R  P  (  p  pop  )  ,  ZRP  (  go  on  )  / I M AG/H I D (  )  ,7  t  p  (  :nr>A] 

COMMON/ DA  T/X  (  -aori  1  j  ,  Y  (  p  r  pp  )  /  T  j  V  /  Pp  R  (  7  .0  n  p  ;  /f.R /F  RrQ  (7  O  pp  )  /fCM^T  /MOPtc 
COw M ON/ S C  -•  L  r  / H *•*  A  X  ,  V ‘  ‘  A  X  / F  I  M r/GJNC 
C 0K'V0\ /  p  OL  /  a  MG  1  3  U  JN  I  T  »  I  G.» MO"  » c*c R  ♦  T 

1 

ION  I  T  =  IJN  I  T 
T  K= I  UN  I T / A 

Al=FLOAT (  I  A  I  )  ‘ 

X  S  C  A  L  =  1 .P/UN  T  T 
!  F  (  K  • c  n • 1  )  X5FAL  =  ^ ./UNIT 
TFK.fq,?  )  v=l  . 

IFU.FQ.l)  V  =  7. 

ZSCAL=V/WAX 

HS  C  A  f_  =  V  /  u  a  x 

IF(<.F0.1>  N=IFIX  f  (rc-p/cfNC)+l  ) 

IFfK'.ro.?)  N  =  NOPTS 

C  ALL  PLOT  (  P  .  p  ,  PV/  dtp  (  1  )  *7. GCAL  *  a  ) 

DO  1  1=1  CJ 

CALL  plot  (FLOAT  (  T  _  1  )  *  X  Sr  a  L  -  °unTP  (  !  )  *?  GC  al  •  ?  ) 

1  COMTINU- 

CALL  GYMpOL ( 0 .n  »? . l  ♦  0 . 1 7  ,  G  H  V  f  R  T  T  C  A  L  .  p  .  p  «  P  ) 

CALL  PLOT (P,n»-P,P«-P ) 

CALL  PLOT ( O.p  «  P  U  n  t  p ( 1  ) *ZGCAL « 0 ) 

DO  7  1  =  1  ,  N 

CALL  PLOT  (FLOAT  (  1-1  )  *X5CAL  iPl.'PTD  (  [  )  *HGCAL  *  ?  ) 

7  CONTI NUF 

CALL  G  v  ‘A  F  0  L  (  n  .  P  C  .  1  C  .  1  ?  1  1  OuHOR  T  70NT  AL  . n, n  .  1 n  j 
Ic ( K.FOc 1 )  CO  TO  A 

CALL  GYM  =  CL  (  FLOAT  (  M!QP  T  f-i  )  ^XGC  AL  +  p  .  G  «  P  .  pp  .  Apcr  r  .  .  p  ,  p  .  ) 

CALL  .GYMRCL  (?  ,P»-1  .G  .0.  IP-^GHTIVf  GYNTHFSIG  Op  TRANSFER  Ft  l.NCT  t  CM  » 

.0.0,15) 

GO  TO  6 

4  DO  G  r=l,M,T< 

CALL  GVMPQL  (  FLOAT  (  j  _  1  )  *  X  GC  A  L  ,  p  .  p  .  p  .  P  p  .  1  p  .  p  .  "  .  -  1  ) 

ArrFRFO  (  I  ) 

CALL  NIP' n  F  P  (  c  L  0  A  T  {  1-1  }  *  X  G  C  AL  ,-O.ig  ,  a  ,  p  .  p  .1  ) 

5  CONTINUE 

CALL  GYM  =  OL  (  FLCAt  (  N  )  *X  GC  \L  +  A  .  G  ,-p  .  P  G  .  1?  .  GH^otO.  ,  p  .  p  ,  g  ) 

IFfL.FQ.!'  CALL  SYMBOL  (  2  . 0  *  -  1  .  S  .  p  ,  i  ?  »  7  RMTPUNC  Atfq  TRawgffp  r(ivrT[C 
.MS *0.0*  ??■  ) 

IF(L.FO.P)  call  SYYHOL ( ? . P » - 1 . G . 1 7 » p 7HNON-TRUNCATFO  TPAMGFFR  F!  !N 

•  C  T  t  onG  ,  o  .  p  C  ?) 

6  COM T  I  NI..IF 

CALL  GYVROL  (?.0^-7.0.0’'.l?»PPMAMGL(r  OF  INCIPFMC^  DFOPFEc,- 

. 0 . 0  .  pp ) 

CALL  N  U  v  p  F  R  (  4 . 1  1-?. D»0, 1  ?  »A  N  G 1.  P»P. o,-l  ) 

CALL  PLOT(0.n,~oc.G,-P) 

WRITE (G*7) 

WRITE (6 *7) 

7  FORMAT ( 1 X 1] PUT  I Mr  IG  OX) 

RETURN 

F  NO 


V 


. 
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SnnpQUTIMF  CUT 

COMMON  /  V  r  P  T  /  D 0  T  p  ('}PO''>)/h0R/P|.'H>TP(  ?  )  /ARFA/qnfl^no) 

rOWMON/nA  T/X  (  poop  )  « Y  f  ^  p  p  0  )  / T  I  v  /  P  c  R  (  7  p  r  p  )  / ^ R  / -  P c Q  (  0  P  o  )  /rON  FT/MOPtc 
COUPON / P  -  a  L  /  HP0  (  Q  ppp  5  »  7  R  p  7  p  p  p  )  /  I Y  a  0  /  h  [  D  (  7  A  p  p  )  ,7{P(7ppp) 
r  o  v m  o  N!  /  f  r  a  l  p  /  h  * '  a  v  ,  v  ‘ 4  *«  x  /  - 1  v  f  /  r.  r  n  c 
rOVMON/RPL/A‘!Glp  »UY  IT  ,  IG*NOp  *tpp  ,T 
on  i  i  =  i*\!nF 
TF  {f  FP--^0(  j  j  j  7,7,3 
7  A=  I  . 

GO  TO  A. 

?  A  =  0. 

4  HP P  {  I  )  =HPP  (  I  1  Tr A 
HTP(  I  )  =  H  T  D  (  j  )-#/*. 

7  R  P (  I ) = 7 P  P ( I) * A 

7  l  p  <  I  )  =  7  T  D  (  I  )  *  A  • 

1  CONTINUE 
WR  I  T  F ( 6 , F  ) 

<=>  FORMAT  (  IX  »oHAFTEP  CUT/} 

P  F  T  u  R 
p  M  0 
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FlJRROUTINF  T CUT 

rOWMON  /  R  PAL  /  HP°  (POOR)  *7RP('*''m)/TMArT/HIP(?ooo) 

COMMON /DAT /XnPn^)-.YnnPO)/TPVPFRf^nno)/rP/FPirOnOPn) /TON  FT /NOPTc 
CON'MON/VErRT/PWnTP(Fcr>0)/HOR/P"nTP('2^oo}/APFA/oiiF(soo) 

COMViON/ROL/  ANGlP  ,llN  I  T  ,  I.G  ,MOc  »CFR  ,  T 
DO  1  1  =  1  ,NDPTF 

TF  (T-ocrp  (f)l  ? 

F  A  =  1  « 

GO  TO  4 
7  A  =  0  • 

4  PV'OTP  (  I  )  =PWDTP  (  I  )  *  A 
PlJOT  P  (  T  }  =  P  u  0  T  P  (  I  )*A 
1  CONTINUF  ■ 

p  r  t  ( i  o  m 
FNin 


4> 


.SUBROUTINE  ^00 ( H * Z * N ) 
REAL  HUH]  ,Z(in),Dno] 

DO  10  I = 1 ,N 

P (  I ) =SORT ( Z (  I ) **?  +  H ( !  )**?) 
10  CONTINUE 
RETURN 
FND 


SUBROUT  I NE  FRFOUM ( F »M , TT  ) 

REAL  FUG) 

00  10  T  =  1,N 

F ( T  )  = FLO AT ( T-l  )  /  ( TT*cLOAT ( M )  ) 
10  CONTINUE 
RFTURN 
END 
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S 1 1 D  P  O  U  T  IMF  TRIJNC(TAI) 


COMMON/ VFRT /PUD TP  (^opn)  /hOR/P"OT  d(  )  /  A  P  E  A  /  D 1 F  (  5  dp  ) 

COMMON /DAT /X  (70  on  )  ,  y  (  ?nPO)  /  t  i  m/PpR  f  3^nn  )  /FR/FRFCf^''^)  /  COM  /NOPt  e 
COMMON  /  RF  AL  /HPP  (  ~  ooo  )  ,^Ro(-5^nA)  /  j  m  a  0  /  H  T  P  (  7  o  n  p  )  ,  7  [  d  /  t  ono  ) 

COMMON  /  S C  A  L  F  /  HM  A  X  ,  V"  A  X  /  -  T  N C  /  0  T  N  r 
COMMON /DOL/ A  MG  IP  »IJN  I  T  *  I  G  *  V0C  ,CFR  ,  t 
C 

IIP  PQPVAT ( 1 X j7 ( 7X » F 1 7 . A ) ) 

116  FORMAT ( ?nx ,p^HREAL  AND  IMAO.  PARTS// ) 

117  PORW AT (  1 X  »  7  ( F 1 0 . 7 , F 1 7 . A  » 1 X  * c 1 7 . A )  ) 

IIS  FORMAT  (  PX  *  7  (  1  X  »  4HTI  vc  ,  *  X  *7HV  — r-OMD.  ,7X  , 7HH-C0MD.  )  )  • 

1 1  Q  FORMAT  (  1  OX  » FHFRcO*  »oy  ,  1AHVPRTTCAL  COMd  .  9  1  7X  ,  1  6HH0P. :  ZONTA  L  COMP. 9 

.  1  7 X  ,  OHV-VODIIH 'S  9  7 X  9  °HH  — MODULUS/  ) 

170  FORMAT ( 1  HI  » 77HTR1 'NCATFD  TR a N ScrR  FUNCTION//) 

171  FORMAT  (  1  HI  9  1  OWT  IMcr  SYNT.//) 

T  7 7  PQRMAT ( 1  HI  9 -HH^AX- 9 r6 . 7 9 IX 9 -HVM AX- . FA. 7 ) 

C 

D  F  L  T  =  1 • /UNIT 
K  K  =  4  *  S  0  0 

CALL  PLOTS  (  oi IF  (  1  )  9^) 

CALL  PLOT ( p . n , 75 . 7 ,-7 ) 

C 

CALL  MAX ( NOP  9  °  U  D  T  D »  H " A  y ) 

CALL  MAX  (  NOP  9  P W DTP  ♦  \/\*  a  X  ) 

WR r  T F ( 6  9 1 7  7 )  HM A X  9 W A  X 
C 

<=1 
L  =  7 

CALL  T I m p (  I  A  I  ,  K 9  L ) 

C 

CALL  CUT 
C 

CALL  ZPPOLNOPTS ,X ) 

CALL  Z F P 0 ( N 0 P T  p  » v ) 

CALL  MOVE ( NOP , 7PD 9 X ) 

CALL  mqvp (NOP 9Z IP  9 Y) 

CALL  LAOS ( NOPTS 9 X ) 

CALL  LAGS(N0PTS,Y1 
N 1 =NOP T  S /? 

CALL  COMCON  (  Y  *  N 1  ) 

CALL  FOUR  m 
CALL  PORT  ID 
CALL  PA C TOO ( X 9 NOPTS ) 

CALL  F A C TOO ( Y  9  N 0 ° T  S ) 

CALL  DtrRICD(°FR  »  NO p  T  S  *  0  I  NC  ) 

CALL  MOVE  (NOPTS  » v  , D D T R  ) 

C 

CALL  ?PRO ( NOPTS 9X ) 

CALL  Z FRO (NOPTS ,v)  • 

CALL  MOVE (NOF ,HRP ,X ) 

CALL  MOVF(NOF,HID9Y) 
call  LAOS (NOPTS  9X) 

CALL  LAOS'f  NOPTS  ,  v  j 

CALL  CO M CON ( Y  9  N 1 ) 


' 
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CALL.  FOUR  ]  D 
CALL  FORT  ID 
CALL  FACTOO ( X 9 NODTS ) 

CALL  p A C TOO ( Y  » NO P t 5 ) 

CALL  MOVF ( MOPTS »X ,PUDTp  )  ( 

C 

c 

WR I T  F ( A *  1 ? 1  ) 

WRITF(6*nR) 

WR I T^ ( A  » 11 7  )  ( PFR (  T  )  »  P W P T P (  I  )  * PUD TP (  I  )  ,  I  =  1  » N l ) 

C 

CALL  MAX (NOPTS*pUDTD ,HMAX) 

CALL  MAX ( MOPTS » PWpTP , VMAX 5 
WR  I  T  E ( 6  *  1 2 2  )  HMAXi V M A  X 
C 

<  =  ? 

L  =  1 

CALL  T ! ME (  I  a  I  ,<,L) 

C 

CALL  T CUT 
C 

CALL  ZFRO ( NOPTS » X ) 

CALL  ZFRO ( NODTS  ,Y ) 

CALL  WOVF  (  MODTS  »PljnTD  ,  X  ) 

C 

CALL  FOUR  ID 
call  FORT  ID 
c 

CAL  L  MOVF ( MOPTS  »  X  »HPD  ) 

C A(L  MOVF (NOPTF  *  V ♦ m I D ) 

CALL  Z  F  d  0 ( N  0  d  T  S  « X ) 

CALL  ZFRO(NODTF»v) 

CALL  MOVc  (MOD  ye  ,  di.iotd  ,  V  )■ 

c 

CALL  FOUR  ID 
CALL  SORT  ID 
C 

CALL  W0VF(\'Odtf*X*ZRP) 

CALL  MO  V  r ( NOD  T  F  »  Y  »  Z  I  °  ) 

CAL.L  MODfZRP  *Z  IP  *  D  WDT  P  »  MOP  T  F  } 

CALL  MOD  {  HR  d  ,HID  ,o!jnjD  ,modtF  ) 

CALL  FRFOUN  (  C‘RFQ  *NOPTF  »  ODLT  ) 

C 

W  R  I  T  F  (  6  ’  1  ?  0  ) 

WR I T F ( 6  » IT  6 ) 

WRITE(6*no) 

WR ITF ( A  » ! ] F )  ( FRFQ (  I  )  ,ZRP ( T )  »Z IP ( T )  >HPD {  I ) ,HJP ( I  >  >dWDTP ( \  )  , 
•  PLJDT  P  (  I  )  *  T  =  1  »NOF  ) 

C 

CALL  max ( NOp  $  PUD TP ,HMAX  ) 

CALL  max  (NOF  ,  d v/ DTd  ,  vmax  ) 

WR  I  TF  (  6  •>  1  ??  )  UMAX,  VmaX 


CALL  T  T M  F ( I  A ! « K  * u 
A  I  =  ( FLOAT ( NOD  T  F ) /UN  I  T ) +? . 
CALL  D!.OT  I  A  I  «  6  •  R  *  -  ^  ) 

R  F  T  UP  N 
FNp 


